math: Remove mpa files (part 2) [BZ #15267]

Previous commit was missing deleted files in sysdeps/ieee754/dbl-64.

Finally remove all mpa related files, headers, declarations, probes, unused
tables and update makefiles.

Reviewed-By: Paul Zimmermann <Paul.Zimmermann@inria.fr>

15 files changed, 0 insertions, 2452 deletions

diff --git a/sysdeps/ieee754/dbl-64/doasin.c b/sysdeps/ieee754/dbl-64/doasin.c
deleted file mode 100644
index a65ef11477..0000000000
--- a/sysdeps/ieee754/dbl-64/doasin.c
+++ /dev/null
@@ -1,81 +0,0 @@
-/*
- * IBM Accurate Mathematical Library
- * written by International Business Machines Corp.
- * Copyright (C) 2001-2021 Free Software Foundation, Inc.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or
- * (at your option) any later version.
- *
- * This program 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program; if not, see <https://www.gnu.org/licenses/>.
- */
-/**********************************************************************/
-/* MODULE_NAME: doasin.c                                              */
-/*                                                                    */
-/* FUNCTION: doasin                                                   */
-/*                                                                    */
-/* FILES NEEDED:endian.h mydefs.h dla.h doasin.h                      */
-/*              mpa.c                                                 */
-/*                                                                    */
-/* Compute arcsin(x,dx,v) of double-length number (x+dx) the result   */
-/* stored in v where v= v[0]+v[1] =arcsin(x+dx)                       */
-/**********************************************************************/
-
-#include "endian.h"
-#include "mydefs.h"
-#include <dla.h>
-#include <math_private.h>
-
-#ifndef SECTION
-# define SECTION
-#endif
-
-/********************************************************************/
-/* Compute arcsin(x,dx,v) of double-length number (x+dx) the result */
-/* stored in v where v= v[0]+v[1] =arcsin(x+dx)                     */
-/********************************************************************/
-void
-SECTION
-__doasin(double x, double dx, double v[]) {
-
-#include "doasin.h"
-
-  static const double
-    d5 =  0.22372159090911789889975459505194491E-01,
-    d6 =  0.17352764422456822913014975683014622E-01,
-    d7 =  0.13964843843786693521653681033981614E-01,
-    d8 =  0.11551791438485242609036067259086589E-01,
-    d9 =  0.97622386568166960207425666787248914E-02,
-    d10 = 0.83638737193775788576092749009744976E-02,
-    d11 = 0.79470250400727425881446981833568758E-02;
-
-  double xx,p,pp,u,uu,r,s;
-  double tc,tcc;
-
-
-/* Taylor series for arcsin for Double-Length numbers         */
-  xx = x*x+2.0*x*dx;
-  p = ((((((d11*xx+d10)*xx+d9)*xx+d8)*xx+d7)*xx+d6)*xx+d5)*xx;
-  pp = 0;
-
-  MUL2(x,dx,x,dx,u,uu,tc,tcc);
-  ADD2(p,pp,c4.x,cc4.x,p,pp,r,s);
-  MUL2(p,pp,u,uu,p,pp,tc,tcc);
-  ADD2(p,pp,c3.x,cc3.x,p,pp,r,s);
-  MUL2(p,pp,u,uu,p,pp,tc,tcc);
-  ADD2(p,pp,c2.x,cc2.x,p,pp,r,s);
-  MUL2(p,pp,u,uu,p,pp,tc,tcc);
-  ADD2(p,pp,c1.x,cc1.x,p,pp,r,s);
-  MUL2(p,pp,u,uu,p,pp,tc,tcc);
-  MUL2(p,pp,x,dx,p,pp,tc,tcc);
-  ADD2(p,pp,x,dx,p,pp,r,s);
-  v[0]=p;
-  v[1]=pp; /* arcsin(x+dx)=v[0]+v[1] */
-}
diff --git a/sysdeps/ieee754/dbl-64/doasin.h b/sysdeps/ieee754/dbl-64/doasin.h
deleted file mode 100644
index bfabca70ae..0000000000
--- a/sysdeps/ieee754/dbl-64/doasin.h
+++ /dev/null
@@ -1,63 +0,0 @@
-
-/*
- * IBM Accurate Mathematical Library
- * Written by International Business Machines Corp.
- * Copyright (C) 2001-2021 Free Software Foundation, Inc.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or
- * (at your option) any later version.
- *
- * This program 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program; if not, see <https://www.gnu.org/licenses/>.
- */
-
-/************************************************************************/
-/*  MODULE_NAME: doasin.h                                                */
-/*                                                                      */
-/*                                                                      */
-/* 	common data and variables definition for BIG or LITTLE ENDIAN   */
-/************************************************************************/
-
-
-
-#ifndef DOASIN_H
-#define DOASIN_H
-
-#ifdef BIG_ENDI
-
- static const  mynumber
-/**/             c1 = {{0x3FC55555, 0x55555555}}, /*  0.16666666666666666    */
-/**/            cc1 = {{0x3C655555, 0x55775389}}, /*  9.2518585419753846e-18 */
-/**/             c2 = {{0x3FB33333, 0x33333333}}, /*  0.074999999999999997   */
-/**/            cc2 = {{0x3C499993, 0x63F1A115}}, /*  2.7755472886508899e-18 */
-/**/             c3 = {{0x3FA6DB6D, 0xB6DB6DB7}}, /*  0.044642857142857144   */
-/**/            cc3 = {{0xBC320FC0, 0x3D5CF0C5}}, /* -9.7911734574147224e-19 */
-/**/             c4 = {{0x3F9F1C71, 0xC71C71C5}}, /*  0.030381944444444437   */
-/**/            cc4 = {{0xBC02B240, 0xFF23ED1E}}; /* -1.2669108566898312e-19 */
-
-#else
-#ifdef LITTLE_ENDI
-
- static const  mynumber
-/**/             c1 = {{0x55555555, 0x3FC55555}}, /*  0.16666666666666666    */
-/**/            cc1 = {{0x55775389, 0x3C655555}}, /*  9.2518585419753846e-18 */
-/**/             c2 = {{0x33333333, 0x3FB33333}}, /*  0.074999999999999997   */
-/**/            cc2 = {{0x63F1A115, 0x3C499993}}, /*  2.7755472886508899e-18 */
-/**/             c3 = {{0xB6DB6DB7, 0x3FA6DB6D}}, /*  0.044642857142857144   */
-/**/            cc3 = {{0x3D5CF0C5, 0xBC320FC0}}, /* -9.7911734574147224e-19 */
-/**/             c4 = {{0xC71C71C5, 0x3F9F1C71}}, /*  0.030381944444444437   */
-/**/            cc4 = {{0xFF23ED1E, 0xBC02B240}}; /* -1.2669108566898312e-19 */
-
-
-#endif
-#endif
-
-
-#endif
diff --git a/sysdeps/ieee754/dbl-64/dosincos.c b/sysdeps/ieee754/dbl-64/dosincos.c
deleted file mode 100644
index 68e3a11401..0000000000
--- a/sysdeps/ieee754/dbl-64/dosincos.c
+++ /dev/null
@@ -1,217 +0,0 @@
-/*
- * IBM Accurate Mathematical Library
- * written by International Business Machines Corp.
- * Copyright (C) 2001-2021 Free Software Foundation, Inc.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or
- * (at your option) any later version.
- *
- * This program 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program; if not, see <https://www.gnu.org/licenses/>.
- */
-/********************************************************************/
-/*                                                                  */
-/* MODULE_NAME: dosincos.c                                          */
-/*                                                                  */
-/*                                                                  */
-/* FUNCTIONS:   dubsin                                              */
-/*              dubcos                                              */
-/*              docos                                               */
-/* FILES NEEDED: endian.h mydefs.h dla.h dosincos.h                 */
-/*               sincos.tbl                                         */
-/*                                                                  */
-/* Routines compute sin() and cos() as Double-Length numbers         */
-/********************************************************************/
-
-
-
-#include "endian.h"
-#include "mydefs.h"
-#include <dla.h>
-#include "dosincos.h"
-#include <math_private.h>
-
-#ifndef SECTION
-# define SECTION
-#endif
-
-extern const union
-{
-  int4 i[880];
-  double x[440];
-} __sincostab attribute_hidden;
-
-/***********************************************************************/
-/* Routine receive Double-Length number (x+dx) and computing sin(x+dx) */
-/* as Double-Length number and store it at array v .It computes it by  */
-/* arithmetic action on Double-Length numbers                          */
-/*(x+dx) between 0 and PI/4                                            */
-/***********************************************************************/
-
-void
-SECTION
-__dubsin (double x, double dx, double v[])
-{
-  double r, s, c, cc, d, dd, d2, dd2, e, ee,
-	 sn, ssn, cs, ccs, ds, dss, dc, dcc;
-  mynumber u;
-  int4 k;
-
-  u.x = x + big.x;
-  k = u.i[LOW_HALF] << 2;
-  x = x - (u.x - big.x);
-  d = x + dx;
-  dd = (x - d) + dx;
-  /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
-  MUL2 (d, dd, d, dd, d2, dd2, c, cc);
-  sn = __sincostab.x[k];       /*                                  */
-  ssn = __sincostab.x[k + 1];  /*      sin(Xi) and cos(Xi)         */
-  cs = __sincostab.x[k + 2];   /*                                  */
-  ccs = __sincostab.x[k + 3];  /*                                  */
-  /* Taylor series for sin ds=sin(t) */
-  MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, c, cc);
-  ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
-  MUL2 (d2, dd2, ds, dss, ds, dss, c, cc);
-  ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
-  MUL2 (d2, dd2, ds, dss, ds, dss, c, cc);
-  MUL2 (d, dd, ds, dss, ds, dss, c, cc);
-  ADD2 (ds, dss, d, dd, ds, dss, r, s);
-
-  /* Taylor series for cos dc=cos(t) */
-  MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, c, cc);
-  ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
-  MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc);
-  ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s);
-  MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc);
-  ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s);
-  MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc);
-
-  MUL2 (cs, ccs, ds, dss, e, ee, c, cc);
-  MUL2 (dc, dcc, sn, ssn, dc, dcc, c, cc);
-  SUB2 (e, ee, dc, dcc, e, ee, r, s);
-  ADD2 (e, ee, sn, ssn, e, ee, r, s);                    /* e+ee=sin(x+dx) */
-
-  v[0] = e;
-  v[1] = ee;
-}
-/**********************************************************************/
-/* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
-/* as Double-Length number and store it in array v .It computes it by */
-/* arithmetic action on Double-Length numbers                         */
-/*(x+dx) between 0 and PI/4                                           */
-/**********************************************************************/
-
-void
-SECTION
-__dubcos (double x, double dx, double v[])
-{
-  double r, s, c, cc, d, dd, d2, dd2, e, ee,
-	 sn, ssn, cs, ccs, ds, dss, dc, dcc;
-  mynumber u;
-  int4 k;
-  u.x = x + big.x;
-  k = u.i[LOW_HALF] << 2;
-  x = x - (u.x - big.x);
-  d = x + dx;
-  dd = (x - d) + dx;  /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
-  MUL2 (d, dd, d, dd, d2, dd2, c, cc);
-  sn = __sincostab.x[k];     /*                                  */
-  ssn = __sincostab.x[k + 1];  /*      sin(Xi) and cos(Xi)         */
-  cs = __sincostab.x[k + 2];   /*                                  */
-  ccs = __sincostab.x[k + 3];  /*                                  */
-  MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, c, cc);
-  ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
-  MUL2 (d2, dd2, ds, dss, ds, dss, c, cc);
-  ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
-  MUL2 (d2, dd2, ds, dss, ds, dss, c, cc);
-  MUL2 (d, dd, ds, dss, ds, dss, c, cc);
-  ADD2 (ds, dss, d, dd, ds, dss, r, s);
-
-  MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, c, cc);
-  ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
-  MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc);
-  ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s);
-  MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc);
-  ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s);
-  MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc);
-
-  MUL2 (cs, ccs, ds, dss, e, ee, c, cc);
-  MUL2 (dc, dcc, sn, ssn, dc, dcc, c, cc);
-
-  MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, c, cc);
-  ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
-  MUL2 (d2, dd2, ds, dss, ds, dss, c, cc);
-  ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
-  MUL2 (d2, dd2, ds, dss, ds, dss, c, cc);
-  MUL2 (d, dd, ds, dss, ds, dss, c, cc);
-  ADD2 (ds, dss, d, dd, ds, dss, r, s);
-  MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, c, cc);
-  ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
-  MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc);
-  ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s);
-  MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc);
-  ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s);
-  MUL2 (d2, dd2, dc, dcc, dc, dcc, c, cc);
-  MUL2 (sn, ssn, ds, dss, e, ee, c, cc);
-  MUL2 (dc, dcc, cs, ccs, dc, dcc, c, cc);
-  ADD2 (e, ee, dc, dcc, e, ee, r, s);
-  SUB2 (cs, ccs, e, ee, e, ee, r, s);
-
-  v[0] = e;
-  v[1] = ee;
-}
-/**********************************************************************/
-/* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
-/* as Double-Length number and store it in array v                    */
-/**********************************************************************/
-void
-SECTION
-__docos (double x, double dx, double v[])
-{
-  double y, yy, p, w[2];
-  if (x > 0)
-    {
-      y = x; yy = dx;
-    }
-  else
-    {
-      y = -x; yy = -dx;
-    }
-  if (y < 0.5 * hp0.x)                                 /*  y< PI/4    */
-    {
-      __dubcos (y, yy, w); v[0] = w[0]; v[1] = w[1];
-    }
-  else if (y < 1.5 * hp0.x)                        /* y< 3/4 * PI */
-    {
-      p = hp0.x - y; /* p = PI/2 - y */
-      yy = hp1.x - yy;
-      y = p + yy;
-      yy = (p - y) + yy;
-      if (y > 0)
-	{
-	  __dubsin (y, yy, w); v[0] = w[0]; v[1] = w[1];
-	}
-      /* cos(x) = sin ( 90 -  x ) */
-      else
-	{
-	  __dubsin (-y, -yy, w); v[0] = -w[0]; v[1] = -w[1];
-	}
-    }
-  else   /* y>= 3/4 * PI */
-    {
-      p = 2.0 * hp0.x - y; /* p = PI- y */
-      yy = 2.0 * hp1.x - yy;
-      y = p + yy;
-      yy = (p - y) + yy;
-      __dubcos (y, yy, w);
-      v[0] = -w[0];
-      v[1] = -w[1];
-    }
-}
diff --git a/sysdeps/ieee754/dbl-64/dosincos.h b/sysdeps/ieee754/dbl-64/dosincos.h
deleted file mode 100644
index 9f34339063..0000000000
--- a/sysdeps/ieee754/dbl-64/dosincos.h
+++ /dev/null
@@ -1,80 +0,0 @@
-
-/*
- * IBM Accurate Mathematical Library
- * Written by International Business Machines Corp.
- * Copyright (C) 2001-2021 Free Software Foundation, Inc.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or
- * (at your option) any later version.
- *
- * This program 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program; if not, see <https://www.gnu.org/licenses/>.
- */
-
-/************************************************************************/
-/*  MODULE_NAME: dosincos.h                                                */
-/*                                                                      */
-/*                                                                      */
-/* 	common data and variables definition for BIG or LITTLE ENDIAN   */
-/************************************************************************/
-
-
-
-#ifndef DOSINCOS_H
-#define DOSINCOS_H
-
-
-#ifdef BIG_ENDI
-static const mynumber
-/**/             s3 = {{0xBFC55555, 0x55555555}},/* -0.16666666666666666    */
-/**/            ss3 = {{0xBC6553AA, 0xE77EE482}},/* -9.2490366677784492e-18 */
-/**/             s5 = {{0x3F811111, 0x11110F15}},/*  0.008333333333332452   */
-/**/            ss5 = {{0xBC21AC06, 0xDA488820}},/* -4.7899996586987931e-19 */
-/**/             s7 = {{0xBF2A019F, 0x5816C78D}},/* -0.00019841261022928957 */
-/**/            ss7 = {{0x3BCDCEC9, 0x6A18BF2A}},/*  1.2624077757871259e-20 */
-/**/             c2 = {{0x3FE00000, 0x00000000}},/*  0.5                    */
-/**/            cc2 = {{0xBA282FD8, 0x00000000}},/* -1.5264073330037701e-28 */
-/**/             c4 = {{0xBFA55555, 0x55555555}},/* -0.041666666666666664   */
-/**/            cc4 = {{0xBC4554BC, 0x2FFF257E}},/* -2.312711276085743e-18  */
-/**/             c6 = {{0x3F56C16C, 0x16C16A96}},/*  0.0013888888888888055  */
-/**/            cc6 = {{0xBBD2E846, 0xE6346F14}},/* -1.6015133010194884e-20 */
-/**/             c8 = {{0xBEFA019F, 0x821D5987}},/* -2.480157866754367e-05  */
-/**/            cc8 = {{0x3B7AB71E, 0x72FFE5CC}},/*  3.5357416224857556e-22 */
-
-/**/            big = {{0x42c80000, 0x00000000}}, /* 52776558133248         */
-
-/**/            hp0 = {{0x3FF921FB, 0x54442D18}}, /* PI / 2                 */
-/**/            hp1 = {{0x3C91A626, 0x33145C07}}; /* 6.123233995736766e-17  */
-#else
-#ifdef LITTLE_ENDI
-static const mynumber
-/**/             s3 = {{0x55555555, 0xBFC55555}},/* -0.16666666666666666    */
-/**/            ss3 = {{0xE77EE482, 0xBC6553AA}},/* -9.2490366677784492e-18 */
-/**/             s5 = {{0x11110F15, 0x3F811111}},/*  0.008333333333332452   */
-/**/            ss5 = {{0xDA488820, 0xBC21AC06}},/* -4.7899996586987931e-19 */
-/**/             s7 = {{0x5816C78D, 0xBF2A019F}},/* -0.00019841261022928957 */
-/**/            ss7 = {{0x6A18BF2A, 0x3BCDCEC9}},/*  1.2624077757871259e-20 */
-/**/             c2 = {{0x00000000, 0x3FE00000}},/*  0.5                    */
-/**/            cc2 = {{0x00000000, 0xBA282FD8}},/* -1.5264073330037701e-28 */
-/**/             c4 = {{0x55555555, 0xBFA55555}},/* -0.041666666666666664   */
-/**/            cc4 = {{0x2FFF257E, 0xBC4554BC}},/* -2.312711276085743e-18  */
-/**/             c6 = {{0x16C16A96, 0x3F56C16C}},/*  0.0013888888888888055  */
-/**/            cc6 = {{0xE6346F14, 0xBBD2E846}},/* -1.6015133010194884e-20 */
-/**/             c8 = {{0x821D5987, 0xBEFA019F}},/* -2.480157866754367e-05  */
-/**/            cc8 = {{0x72FFE5CC, 0x3B7AB71E}},/*  3.5357416224857556e-22 */
-
-/**/            big = {{0x00000000, 0x42c80000}}, /* 52776558133248         */
-
-/**/            hp0 = {{0x54442D18, 0x3FF921FB}}, /* PI / 2                 */
-/**/            hp1 = {{0x33145C07, 0x3C91A626}}; /* 6.123233995736766e-17  */
-#endif
-#endif
-
-#endif
diff --git a/sysdeps/ieee754/dbl-64/mpa-arch.h b/sysdeps/ieee754/dbl-64/mpa-arch.h
deleted file mode 100644
index fbe296d8b5..0000000000
--- a/sysdeps/ieee754/dbl-64/mpa-arch.h
+++ /dev/null
@@ -1,47 +0,0 @@
-/* Overridable constants and operations.
-   Copyright (C) 2013-2021 Free Software Foundation, Inc.
-
-   This program is free software; you can redistribute it and/or modify
-   it under the terms of the GNU Lesser General Public License as published by
-   the Free Software Foundation; either version 2.1 of the License, or
-   (at your option) any later version.
-
-   This program 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 Lesser General Public License for more details.
-
-   You should have received a copy of the GNU Lesser General Public License
-   along with this program; if not, see <https://www.gnu.org/licenses/>.  */
-
-#include <stdint.h>
-
-typedef long mantissa_t;
-typedef int64_t mantissa_store_t;
-
-#define TWOPOW(i) (1L << i)
-
-#define RADIX_EXP 24
-#define  RADIX TWOPOW (RADIX_EXP)               /* 2^24    */
-
-/* Divide D by RADIX and put the remainder in R.  D must be a non-negative
-   integral value.  */
-#define DIV_RADIX(d, r) \
-  ({                                                                         \
-     r = d & (RADIX - 1);                                                    \
-     d >>= RADIX_EXP;                                                        \
-   })
-
-/* Put the integer component of a double X in R and retain the fraction in
-   X.  This is used in extracting mantissa digits for MP_NO by using the
-   integer portion of the current value of the number as the current mantissa
-   digit and then scaling by RADIX to get the next mantissa digit in the same
-   manner.  */
-#define INTEGER_OF(x, i) \
-  ({                                                                          \
-     i = (mantissa_t) x;                                                       \
-     x -= i;                                                                   \
-   })
-
-/* Align IN down to F.  The code assumes that F is a power of two.  */
-#define ALIGN_DOWN_TO(in, f) ((in) & - (f))
diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c
deleted file mode 100644
index eb5d8e8e89..0000000000
--- a/sysdeps/ieee754/dbl-64/mpa.c
+++ /dev/null
@@ -1,913 +0,0 @@
-/*
- * IBM Accurate Mathematical Library
- * written by International Business Machines Corp.
- * Copyright (C) 2001-2021 Free Software Foundation, Inc.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or
- * (at your option) any later version.
- *
- * This program 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program; if not, see <https://www.gnu.org/licenses/>.
- */
-/************************************************************************/
-/*  MODULE_NAME: mpa.c                                                  */
-/*                                                                      */
-/*  FUNCTIONS:                                                          */
-/*               mcr                                                    */
-/*               acr                                                    */
-/*               cpy                                                    */
-/*               norm                                                   */
-/*               denorm                                                 */
-/*               mp_dbl                                                 */
-/*               dbl_mp                                                 */
-/*               add_magnitudes                                         */
-/*               sub_magnitudes                                         */
-/*               add                                                    */
-/*               sub                                                    */
-/*               mul                                                    */
-/*               inv                                                    */
-/*               dvd                                                    */
-/*                                                                      */
-/* Arithmetic functions for multiple precision numbers.                 */
-/* Relative errors are bounded                                          */
-/************************************************************************/
-
-
-#include "endian.h"
-#include "mpa.h"
-#include <sys/param.h>
-#include <alloca.h>
-
-#ifndef SECTION
-# define SECTION
-#endif
-
-#ifndef NO__CONST
-const mp_no __mpone = { 1, { 1.0, 1.0 } };
-const mp_no __mptwo = { 1, { 1.0, 2.0 } };
-#endif
-
-#ifndef NO___ACR
-/* Compare mantissa of two multiple precision numbers regardless of the sign
-   and exponent of the numbers.  */
-static int
-mcr (const mp_no *x, const mp_no *y, int p)
-{
-  long i;
-  long p2 = p;
-  for (i = 1; i <= p2; i++)
-    {
-      if (X[i] == Y[i])
-	continue;
-      else if (X[i] > Y[i])
-	return 1;
-      else
-	return -1;
-    }
-  return 0;
-}
-
-/* Compare the absolute values of two multiple precision numbers.  */
-int
-__acr (const mp_no *x, const mp_no *y, int p)
-{
-  long i;
-
-  if (X[0] == 0)
-    {
-      if (Y[0] == 0)
-	i = 0;
-      else
-	i = -1;
-    }
-  else if (Y[0] == 0)
-    i = 1;
-  else
-    {
-      if (EX > EY)
-	i = 1;
-      else if (EX < EY)
-	i = -1;
-      else
-	i = mcr (x, y, p);
-    }
-
-  return i;
-}
-#endif
-
-#ifndef NO___CPY
-/* Copy multiple precision number X into Y.  They could be the same
-   number.  */
-void
-__cpy (const mp_no *x, mp_no *y, int p)
-{
-  long i;
-
-  EY = EX;
-  for (i = 0; i <= p; i++)
-    Y[i] = X[i];
-}
-#endif
-
-#ifndef NO___MP_DBL
-/* Convert a multiple precision number *X into a double precision
-   number *Y, normalized case (|x| >= 2**(-1022))).  X has precision
-   P, which is positive.  */
-static void
-norm (const mp_no *x, double *y, int p)
-{
-# define R RADIXI
-  long i;
-  double c;
-  mantissa_t a, u, v, z[5];
-  if (p < 5)
-    {
-      if (p == 1)
-	c = X[1];
-      else if (p == 2)
-	c = X[1] + R * X[2];
-      else if (p == 3)
-	c = X[1] + R * (X[2] + R * X[3]);
-      else /* p == 4.  */
-	c = (X[1] + R * X[2]) + R * R * (X[3] + R * X[4]);
-    }
-  else
-    {
-      for (a = 1, z[1] = X[1]; z[1] < TWO23; )
-	{
-	  a *= 2;
-	  z[1] *= 2;
-	}
-
-      for (i = 2; i < 5; i++)
-	{
-	  mantissa_store_t d, r;
-	  d = X[i] * (mantissa_store_t) a;
-	  DIV_RADIX (d, r);
-	  z[i] = r;
-	  z[i - 1] += d;
-	}
-
-      u = ALIGN_DOWN_TO (z[3], TWO19);
-      v = z[3] - u;
-
-      if (v == TWO18)
-	{
-	  if (z[4] == 0)
-	    {
-	      for (i = 5; i <= p; i++)
-		{
-		  if (X[i] == 0)
-		    continue;
-		  else
-		    {
-		      z[3] += 1;
-		      break;
-		    }
-		}
-	    }
-	  else
-	    z[3] += 1;
-	}
-
-      c = (z[1] + R * (z[2] + R * z[3])) / a;
-    }
-
-  c *= X[0];
-
-  for (i = 1; i < EX; i++)
-    c *= RADIX;
-  for (i = 1; i > EX; i--)
-    c *= RADIXI;
-
-  *y = c;
-# undef R
-}
-
-/* Convert a multiple precision number *X into a double precision
-   number *Y, Denormal case  (|x| < 2**(-1022))).  */
-static void
-denorm (const mp_no *x, double *y, int p)
-{
-  long i, k;
-  long p2 = p;
-  double c;
-  mantissa_t u, z[5];
-
-# define R RADIXI
-  if (EX < -44 || (EX == -44 && X[1] < TWO5))
-    {
-      *y = 0;
-      return;
-    }
-
-  if (p2 == 1)
-    {
-      if (EX == -42)
-	{
-	  z[1] = X[1] + TWO10;
-	  z[2] = 0;
-	  z[3] = 0;
-	  k = 3;
-	}
-      else if (EX == -43)
-	{
-	  z[1] = TWO10;
-	  z[2] = X[1];
-	  z[3] = 0;
-	  k = 2;
-	}
-      else
-	{
-	  z[1] = TWO10;
-	  z[2] = 0;
-	  z[3] = X[1];
-	  k = 1;
-	}
-    }
-  else if (p2 == 2)
-    {
-      if (EX == -42)
-	{
-	  z[1] = X[1] + TWO10;
-	  z[2] = X[2];
-	  z[3] = 0;
-	  k = 3;
-	}
-      else if (EX == -43)
-	{
-	  z[1] = TWO10;
-	  z[2] = X[1];
-	  z[3] = X[2];
-	  k = 2;
-	}
-      else
-	{
-	  z[1] = TWO10;
-	  z[2] = 0;
-	  z[3] = X[1];
-	  k = 1;
-	}
-    }
-  else
-    {
-      if (EX == -42)
-	{
-	  z[1] = X[1] + TWO10;
-	  z[2] = X[2];
-	  k = 3;
-	}
-      else if (EX == -43)
-	{
-	  z[1] = TWO10;
-	  z[2] = X[1];
-	  k = 2;
-	}
-      else
-	{
-	  z[1] = TWO10;
-	  z[2] = 0;
-	  k = 1;
-	}
-      z[3] = X[k];
-    }
-
-  u = ALIGN_DOWN_TO (z[3], TWO5);
-
-  if (u == z[3])
-    {
-      for (i = k + 1; i <= p2; i++)
-	{
-	  if (X[i] == 0)
-	    continue;
-	  else
-	    {
-	      z[3] += 1;
-	      break;
-	    }
-	}
-    }
-
-  c = X[0] * ((z[1] + R * (z[2] + R * z[3])) - TWO10);
-
-  *y = c * TWOM1032;
-# undef R
-}
-
-/* Convert multiple precision number *X into double precision number *Y.  The
-   result is correctly rounded to the nearest/even.  */
-void
-__mp_dbl (const mp_no *x, double *y, int p)
-{
-  if (X[0] == 0)
-    {
-      *y = 0;
-      return;
-    }
-
-  if (__glibc_likely (EX > -42 || (EX == -42 && X[1] >= TWO10)))
-    norm (x, y, p);
-  else
-    denorm (x, y, p);
-}
-#endif
-
-/* Get the multiple precision equivalent of X into *Y.  If the precision is too
-   small, the result is truncated.  */
-void
-SECTION
-__dbl_mp (double x, mp_no *y, int p)
-{
-  long i, n;
-  long p2 = p;
-
-  /* Sign.  */
-  if (x == 0)
-    {
-      Y[0] = 0;
-      return;
-    }
-  else if (x > 0)
-    Y[0] = 1;
-  else
-    {
-      Y[0] = -1;
-      x = -x;
-    }
-
-  /* Exponent.  */
-  for (EY = 1; x >= RADIX; EY += 1)
-    x *= RADIXI;
-  for (; x < 1; EY -= 1)
-    x *= RADIX;
-
-  /* Digits.  */
-  n = MIN (p2, 4);
-  for (i = 1; i <= n; i++)
-    {
-      INTEGER_OF (x, Y[i]);
-      x *= RADIX;
-    }
-  for (; i <= p2; i++)
-    Y[i] = 0;
-}
-
-/* Add magnitudes of *X and *Y assuming that abs (*X) >= abs (*Y) > 0.  The
-   sign of the sum *Z is not changed.  X and Y may overlap but not X and Z or
-   Y and Z.  No guard digit is used.  The result equals the exact sum,
-   truncated.  */
-static void
-SECTION
-add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  long i, j, k;
-  long p2 = p;
-  mantissa_t zk;
-
-  EZ = EX;
-
-  i = p2;
-  j = p2 + EY - EX;
-  k = p2 + 1;
-
-  if (__glibc_unlikely (j < 1))
-    {
-      __cpy (x, z, p);
-      return;
-    }
-
-  zk = 0;
-
-  for (; j > 0; i--, j--)
-    {
-      zk += X[i] + Y[j];
-      if (zk >= RADIX)
-	{
-	  Z[k--] = zk - RADIX;
-	  zk = 1;
-	}
-      else
-	{
-	  Z[k--] = zk;
-	  zk = 0;
-	}
-    }
-
-  for (; i > 0; i--)
-    {
-      zk += X[i];
-      if (zk >= RADIX)
-	{
-	  Z[k--] = zk - RADIX;
-	  zk = 1;
-	}
-      else
-	{
-	  Z[k--] = zk;
-	  zk = 0;
-	}
-    }
-
-  if (zk == 0)
-    {
-      for (i = 1; i <= p2; i++)
-	Z[i] = Z[i + 1];
-    }
-  else
-    {
-      Z[1] = zk;
-      EZ += 1;
-    }
-}
-
-/* Subtract the magnitudes of *X and *Y assuming that abs (*x) > abs (*y) > 0.
-   The sign of the difference *Z is not changed.  X and Y may overlap but not X
-   and Z or Y and Z.  One guard digit is used.  The error is less than one
-   ULP.  */
-static void
-SECTION
-sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  long i, j, k;
-  long p2 = p;
-  mantissa_t zk;
-
-  EZ = EX;
-  i = p2;
-  j = p2 + EY - EX;
-  k = p2;
-
-  /* Y is too small compared to X, copy X over to the result.  */
-  if (__glibc_unlikely (j < 1))
-    {
-      __cpy (x, z, p);
-      return;
-    }
-
-  /* The relevant least significant digit in Y is non-zero, so we factor it in
-     to enhance accuracy.  */
-  if (j < p2 && Y[j + 1] > 0)
-    {
-      Z[k + 1] = RADIX - Y[j + 1];
-      zk = -1;
-    }
-  else
-    zk = Z[k + 1] = 0;
-
-  /* Subtract and borrow.  */
-  for (; j > 0; i--, j--)
-    {
-      zk += (X[i] - Y[j]);
-      if (zk < 0)
-	{
-	  Z[k--] = zk + RADIX;
-	  zk = -1;
-	}
-      else
-	{
-	  Z[k--] = zk;
-	  zk = 0;
-	}
-    }
-
-  /* We're done with digits from Y, so it's just digits in X.  */
-  for (; i > 0; i--)
-    {
-      zk += X[i];
-      if (zk < 0)
-	{
-	  Z[k--] = zk + RADIX;
-	  zk = -1;
-	}
-      else
-	{
-	  Z[k--] = zk;
-	  zk = 0;
-	}
-    }
-
-  /* Normalize.  */
-  for (i = 1; Z[i] == 0; i++)
-    ;
-  EZ = EZ - i + 1;
-  for (k = 1; i <= p2 + 1; )
-    Z[k++] = Z[i++];
-  for (; k <= p2; )
-    Z[k++] = 0;
-}
-
-/* Add *X and *Y and store the result in *Z.  X and Y may overlap, but not X
-   and Z or Y and Z.  One guard digit is used.  The error is less than one
-   ULP.  */
-void
-SECTION
-__add (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  int n;
-
-  if (X[0] == 0)
-    {
-      __cpy (y, z, p);
-      return;
-    }
-  else if (Y[0] == 0)
-    {
-      __cpy (x, z, p);
-      return;
-    }
-
-  if (X[0] == Y[0])
-    {
-      if (__acr (x, y, p) > 0)
-	{
-	  add_magnitudes (x, y, z, p);
-	  Z[0] = X[0];
-	}
-      else
-	{
-	  add_magnitudes (y, x, z, p);
-	  Z[0] = Y[0];
-	}
-    }
-  else
-    {
-      if ((n = __acr (x, y, p)) == 1)
-	{
-	  sub_magnitudes (x, y, z, p);
-	  Z[0] = X[0];
-	}
-      else if (n == -1)
-	{
-	  sub_magnitudes (y, x, z, p);
-	  Z[0] = Y[0];
-	}
-      else
-	Z[0] = 0;
-    }
-}
-
-/* Subtract *Y from *X and return the result in *Z.  X and Y may overlap but
-   not X and Z or Y and Z.  One guard digit is used.  The error is less than
-   one ULP.  */
-void
-SECTION
-__sub (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  int n;
-
-  if (X[0] == 0)
-    {
-      __cpy (y, z, p);
-      Z[0] = -Z[0];
-      return;
-    }
-  else if (Y[0] == 0)
-    {
-      __cpy (x, z, p);
-      return;
-    }
-
-  if (X[0] != Y[0])
-    {
-      if (__acr (x, y, p) > 0)
-	{
-	  add_magnitudes (x, y, z, p);
-	  Z[0] = X[0];
-	}
-      else
-	{
-	  add_magnitudes (y, x, z, p);
-	  Z[0] = -Y[0];
-	}
-    }
-  else
-    {
-      if ((n = __acr (x, y, p)) == 1)
-	{
-	  sub_magnitudes (x, y, z, p);
-	  Z[0] = X[0];
-	}
-      else if (n == -1)
-	{
-	  sub_magnitudes (y, x, z, p);
-	  Z[0] = -Y[0];
-	}
-      else
-	Z[0] = 0;
-    }
-}
-
-#ifndef NO__MUL
-/* Multiply *X and *Y and store result in *Z.  X and Y may overlap but not X
-   and Z or Y and Z.  For P in [1, 2, 3], the exact result is truncated to P
-   digits.  In case P > 3 the error is bounded by 1.001 ULP.  */
-void
-SECTION
-__mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  long i, j, k, ip, ip2;
-  long p2 = p;
-  mantissa_store_t zk;
-  const mp_no *a;
-  mantissa_store_t *diag;
-
-  /* Is z=0?  */
-  if (__glibc_unlikely (X[0] * Y[0] == 0))
-    {
-      Z[0] = 0;
-      return;
-    }
-
-  /* We need not iterate through all X's and Y's since it's pointless to
-     multiply zeroes.  Here, both are zero...  */
-  for (ip2 = p2; ip2 > 0; ip2--)
-    if (X[ip2] != 0 || Y[ip2] != 0)
-      break;
-
-  a = X[ip2] != 0 ? y : x;
-
-  /* ... and here, at least one of them is still zero.  */
-  for (ip = ip2; ip > 0; ip--)
-    if (a->d[ip] != 0)
-      break;
-
-  /* The product looks like this for p = 3 (as an example):
-
-
-				a1    a2    a3
-		 x		b1    b2    b3
-		 -----------------------------
-			     a1*b3 a2*b3 a3*b3
-		       a1*b2 a2*b2 a3*b2
-		 a1*b1 a2*b1 a3*b1
-
-     So our K needs to ideally be P*2, but we're limiting ourselves to P + 3
-     for P >= 3.  We compute the above digits in two parts; the last P-1
-     digits and then the first P digits.  The last P-1 digits are a sum of
-     products of the input digits from P to P-k where K is 0 for the least
-     significant digit and increases as we go towards the left.  The product
-     term is of the form X[k]*X[P-k] as can be seen in the above example.
-
-     The first P digits are also a sum of products with the same product term,
-     except that the sum is from 1 to k.  This is also evident from the above
-     example.
-
-     Another thing that becomes evident is that only the most significant
-     ip+ip2 digits of the result are non-zero, where ip and ip2 are the
-     'internal precision' of the input numbers, i.e. digits after ip and ip2
-     are all 0.  */
-
-  k = (__glibc_unlikely (p2 < 3)) ? p2 + p2 : p2 + 3;
-
-  while (k > ip + ip2 + 1)
-    Z[k--] = 0;
-
-  zk = 0;
-
-  /* Precompute sums of diagonal elements so that we can directly use them
-     later.  See the next comment to know we why need them.  */
-  diag = alloca (k * sizeof (mantissa_store_t));
-  mantissa_store_t d = 0;
-  for (i = 1; i <= ip; i++)
-    {
-      d += X[i] * (mantissa_store_t) Y[i];
-      diag[i] = d;
-    }
-  while (i < k)
-    diag[i++] = d;
-
-  while (k > p2)
-    {
-      long lim = k / 2;
-
-      if (k % 2 == 0)
-	/* We want to add this only once, but since we subtract it in the sum
-	   of products above, we add twice.  */
-	zk += 2 * X[lim] * (mantissa_store_t) Y[lim];
-
-      for (i = k - p2, j = p2; i < j; i++, j--)
-	zk += (X[i] + X[j]) * (mantissa_store_t) (Y[i] + Y[j]);
-
-      zk -= diag[k - 1];
-
-      DIV_RADIX (zk, Z[k]);
-      k--;
-    }
-
-  /* The real deal.  Mantissa digit Z[k] is the sum of all X[i] * Y[j] where i
-     goes from 1 -> k - 1 and j goes the same range in reverse.  To reduce the
-     number of multiplications, we halve the range and if k is an even number,
-     add the diagonal element X[k/2]Y[k/2].  Through the half range, we compute
-     X[i] * Y[j] as (X[i] + X[j]) * (Y[i] + Y[j]) - X[i] * Y[i] - X[j] * Y[j].
-
-     This reduction tells us that we're summing two things, the first term
-     through the half range and the negative of the sum of the product of all
-     terms of X and Y in the full range.  i.e.
-
-     SUM(X[i] * Y[i]) for k terms.  This is precalculated above for each k in
-     a single loop so that it completes in O(n) time and can hence be directly
-     used in the loop below.  */
-  while (k > 1)
-    {
-      long lim = k / 2;
-
-      if (k % 2 == 0)
-	/* We want to add this only once, but since we subtract it in the sum
-	   of products above, we add twice.  */
-        zk += 2 * X[lim] * (mantissa_store_t) Y[lim];
-
-      for (i = 1, j = k - 1; i < j; i++, j--)
-	zk += (X[i] + X[j]) * (mantissa_store_t) (Y[i] + Y[j]);
-
-      zk -= diag[k - 1];
-
-      DIV_RADIX (zk, Z[k]);
-      k--;
-    }
-  Z[k] = zk;
-
-  /* Get the exponent sum into an intermediate variable.  This is a subtle
-     optimization, where given enough registers, all operations on the exponent
-     happen in registers and the result is written out only once into EZ.  */
-  int e = EX + EY;
-
-  /* Is there a carry beyond the most significant digit?  */
-  if (__glibc_unlikely (Z[1] == 0))
-    {
-      for (i = 1; i <= p2; i++)
-	Z[i] = Z[i + 1];
-      e--;
-    }
-
-  EZ = e;
-  Z[0] = X[0] * Y[0];
-}
-#endif
-
-#ifndef NO__SQR
-/* Square *X and store result in *Y.  X and Y may not overlap.  For P in
-   [1, 2, 3], the exact result is truncated to P digits.  In case P > 3 the
-   error is bounded by 1.001 ULP.  This is a faster special case of
-   multiplication.  */
-void
-SECTION
-__sqr (const mp_no *x, mp_no *y, int p)
-{
-  long i, j, k, ip;
-  mantissa_store_t yk;
-
-  /* Is z=0?  */
-  if (__glibc_unlikely (X[0] == 0))
-    {
-      Y[0] = 0;
-      return;
-    }
-
-  /* We need not iterate through all X's since it's pointless to
-     multiply zeroes.  */
-  for (ip = p; ip > 0; ip--)
-    if (X[ip] != 0)
-      break;
-
-  k = (__glibc_unlikely (p < 3)) ? p + p : p + 3;
-
-  while (k > 2 * ip + 1)
-    Y[k--] = 0;
-
-  yk = 0;
-
-  while (k > p)
-    {
-      mantissa_store_t yk2 = 0;
-      long lim = k / 2;
-
-      if (k % 2 == 0)
-	yk += X[lim] * (mantissa_store_t) X[lim];
-
-      /* In __mul, this loop (and the one within the next while loop) run
-         between a range to calculate the mantissa as follows:
-
-         Z[k] = X[k] * Y[n] + X[k+1] * Y[n-1] ... + X[n-1] * Y[k+1]
-		+ X[n] * Y[k]
-
-         For X == Y, we can get away with summing halfway and doubling the
-	 result.  For cases where the range size is even, the mid-point needs
-	 to be added separately (above).  */
-      for (i = k - p, j = p; i < j; i++, j--)
-	yk2 += X[i] * (mantissa_store_t) X[j];
-
-      yk += 2 * yk2;
-
-      DIV_RADIX (yk, Y[k]);
-      k--;
-    }
-
-  while (k > 1)
-    {
-      mantissa_store_t yk2 = 0;
-      long lim = k / 2;
-
-      if (k % 2 == 0)
-	yk += X[lim] * (mantissa_store_t) X[lim];
-
-      /* Likewise for this loop.  */
-      for (i = 1, j = k - 1; i < j; i++, j--)
-	yk2 += X[i] * (mantissa_store_t) X[j];
-
-      yk += 2 * yk2;
-
-      DIV_RADIX (yk, Y[k]);
-      k--;
-    }
-  Y[k] = yk;
-
-  /* Squares are always positive.  */
-  Y[0] = 1;
-
-  /* Get the exponent sum into an intermediate variable.  This is a subtle
-     optimization, where given enough registers, all operations on the exponent
-     happen in registers and the result is written out only once into EZ.  */
-  int e = EX * 2;
-
-  /* Is there a carry beyond the most significant digit?  */
-  if (__glibc_unlikely (Y[1] == 0))
-    {
-      for (i = 1; i <= p; i++)
-	Y[i] = Y[i + 1];
-      e--;
-    }
-
-  EY = e;
-}
-#endif
-
-/* Invert *X and store in *Y.  Relative error bound:
-   - For P = 2: 1.001 * R ^ (1 - P)
-   - For P = 3: 1.063 * R ^ (1 - P)
-   - For P > 3: 2.001 * R ^ (1 - P)
-
-   *X = 0 is not permissible.  */
-static void
-SECTION
-__inv (const mp_no *x, mp_no *y, int p)
-{
-  long i;
-  double t;
-  mp_no z, w;
-  static const int np1[] =
-    { 0, 0, 0, 0, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3,
-    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4
-  };
-
-  __cpy (x, &z, p);
-  z.e = 0;
-  __mp_dbl (&z, &t, p);
-  t = 1 / t;
-
-  /* t == 0 will never happen at this point, since 1/t can only be 0 if t is
-     infinity, but before the division t == mantissa of x (exponent is 0).  We
-     can instruct the compiler to ignore this case.  */
-  if (t == 0)
-    __builtin_unreachable ();
-
-  __dbl_mp (t, y, p);
-  EY -= EX;
-
-  for (i = 0; i < np1[p]; i++)
-    {
-      __cpy (y, &w, p);
-      __mul (x, &w, y, p);
-      __sub (&__mptwo, y, &z, p);
-      __mul (&w, &z, y, p);
-    }
-}
-
-/* Divide *X by *Y and store result in *Z.  X and Y may overlap but not X and Z
-   or Y and Z.  Relative error bound:
-   - For P = 2: 2.001 * R ^ (1 - P)
-   - For P = 3: 2.063 * R ^ (1 - P)
-   - For P > 3: 3.001 * R ^ (1 - P)
-
-   *X = 0 is not permissible.  */
-void
-SECTION
-__dvd (const mp_no *x, const mp_no *y, mp_no *z, int p)
-{
-  mp_no w;
-
-  if (X[0] == 0)
-    Z[0] = 0;
-  else
-    {
-      __inv (y, &w, p);
-      __mul (x, &w, z, p);
-    }
-}
diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h
deleted file mode 100644
index c28630e148..0000000000
--- a/sysdeps/ieee754/dbl-64/mpa.h
+++ /dev/null
@@ -1,123 +0,0 @@
-/*
- * IBM Accurate Mathematical Library
- * Written by International Business Machines Corp.
- * Copyright (C) 2001-2021 Free Software Foundation, Inc.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or
- * (at your option) any later version.
- *
- * This program 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program; if not, see <https://www.gnu.org/licenses/>.
- */
-
-/************************************************************************/
-/*  MODULE_NAME: mpa.h                                                  */
-/*                                                                      */
-/*  FUNCTIONS:                                                          */
-/*               mcr                                                    */
-/*               acr                                                    */
-/*               cpy                                                    */
-/*               mp_dbl                                                 */
-/*               dbl_mp                                                 */
-/*               add                                                    */
-/*               sub                                                    */
-/*               mul                                                    */
-/*               dvd                                                    */
-/*                                                                      */
-/* Arithmetic functions for multiple precision numbers.                 */
-/* Common types and definition                                          */
-/************************************************************************/
-
-#include <mpa-arch.h>
-
-/* The mp_no structure holds the details of a multi-precision floating point
-   number.
-
-   - The radix of the number (R) is 2 ^ 24.
-
-   - E: The exponent of the number.
-
-   - D[0]: The sign (-1, 1) or 0 if the value is 0.  In the latter case, the
-     values of the remaining members of the structure are ignored.
-
-   - D[1] - D[p]: The mantissa of the number where:
-
-	0 <= D[i] < R and
-	P is the precision of the number and 1 <= p <= 32
-
-     D[p+1] ... D[39] have no significance.
-
-   - The value of the number is:
-
-	D[1] * R ^ (E - 1) + D[2] * R ^ (E - 2) ... D[p] * R ^ (E - p)
-
-   */
-typedef struct
-{
-  int e;
-  mantissa_t d[40];
-} mp_no;
-
-typedef union
-{
-  int i[2];
-  double d;
-} number;
-
-extern const mp_no __mpone;
-extern const mp_no __mptwo;
-
-#define  X   x->d
-#define  Y   y->d
-#define  Z   z->d
-#define  EX  x->e
-#define  EY  y->e
-#define  EZ  z->e
-
-#ifndef RADIXI
-# define  RADIXI    0x1.0p-24		/* 2^-24   */
-#endif
-
-#ifndef TWO52
-# define  TWO52     0x1.0p52		/* 2^52    */
-#endif
-
-#define  TWO5      TWOPOW (5)		/* 2^5     */
-#define  TWO8      TWOPOW (8)		/* 2^52    */
-#define  TWO10     TWOPOW (10)		/* 2^10    */
-#define  TWO18     TWOPOW (18)		/* 2^18    */
-#define  TWO19     TWOPOW (19)		/* 2^19    */
-#define  TWO23     TWOPOW (23)		/* 2^23    */
-
-#define  HALFRAD   TWO23
-
-#define  TWO57     0x1.0p57		/* 2^57    */
-#define  TWO71     0x1.0p71		/* 2^71    */
-#define  TWOM1032  0x1.0p-1032		/* 2^-1032 */
-#define  TWOM1022  0x1.0p-1022		/* 2^-1022 */
-
-#define  HALF      0x1.0p-1		/* 1/2 */
-#define  MHALF     -0x1.0p-1		/* -1/2 */
-
-int __acr (const mp_no *, const mp_no *, int);
-void __cpy (const mp_no *, mp_no *, int);
-void __mp_dbl (const mp_no *, double *, int);
-void __dbl_mp (double, mp_no *, int);
-void __add (const mp_no *, const mp_no *, mp_no *, int);
-void __sub (const mp_no *, const mp_no *, mp_no *, int);
-void __mul (const mp_no *, const mp_no *, mp_no *, int);
-void __sqr (const mp_no *, mp_no *, int);
-void __dvd (const mp_no *, const mp_no *, mp_no *, int);
-
-extern void __mpatan (mp_no *, mp_no *, int);
-extern void __mpatan2 (mp_no *, mp_no *, mp_no *, int);
-extern void __mpsqrt (mp_no *, mp_no *, int);
-extern void __c32 (mp_no *, mp_no *, mp_no *, int);
-extern int __mpranred (double, mp_no *, int);
diff --git a/sysdeps/ieee754/dbl-64/mpatan.c b/sysdeps/ieee754/dbl-64/mpatan.c
deleted file mode 100644
index 77f84ddbf2..0000000000
--- a/sysdeps/ieee754/dbl-64/mpatan.c
+++ /dev/null
@@ -1,116 +0,0 @@
-/*
- * IBM Accurate Mathematical Library
- * written by International Business Machines Corp.
- * Copyright (C) 2001-2021 Free Software Foundation, Inc.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or
- * (at your option) any later version.
- *
- * This program 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU  Lesser General Public License
- * along with this program; if not, see <https://www.gnu.org/licenses/>.
- */
-/******************************************************************/
-/*                                                                */
-/* MODULE_NAME:mpatan.c                                           */
-/*                                                                */
-/* FUNCTIONS:mpatan                                               */
-/*                                                                */
-/* FILES NEEDED: mpa.h endian.h mpatan.h                          */
-/*               mpa.c                                            */
-/*                                                                */
-/* Multi-Precision Atan function subroutine, for precision p >= 4.*/
-/* The relative error of the result is bounded by 34.32*r**(1-p), */
-/* where r=2**24.                                                 */
-/******************************************************************/
-
-#include "endian.h"
-#include "mpa.h"
-#include <math.h>
-
-#ifndef SECTION
-# define SECTION
-#endif
-
-#include "mpatan.h"
-
-void
-SECTION
-__mpatan (mp_no *x, mp_no *y, int p)
-{
-  int i, m, n;
-  double dx;
-  mp_no mptwoim1 =
-  {
-    0,
-    {
-      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
-    }
-  };
-
-  mp_no mps, mpsm, mpt, mpt1, mpt2, mpt3;
-
-  /* Choose m and initiate mptwoim1.  */
-  if (EX > 0)
-    m = 7;
-  else if (EX < 0)
-    m = 0;
-  else
-    {
-      __mp_dbl (x, &dx, p);
-      dx = fabs (dx);
-      for (m = 6; m > 0; m--)
-	{
-	  if (dx > __atan_xm[m].d)
-	    break;
-	}
-    }
-  mptwoim1.e = 1;
-  mptwoim1.d[0] = 1;
-
-  /* Reduce x m times.  */
-  __sqr (x, &mpsm, p);
-  if (m == 0)
-    __cpy (x, &mps, p);
-  else
-    {
-      for (i = 0; i < m; i++)
-	{
-	  __add (&__mpone, &mpsm, &mpt1, p);
-	  __mpsqrt (&mpt1, &mpt2, p);
-	  __add (&mpt2, &mpt2, &mpt1, p);
-	  __add (&__mptwo, &mpsm, &mpt2, p);
-	  __add (&mpt1, &mpt2, &mpt3, p);
-	  __dvd (&mpsm, &mpt3, &mpt1, p);
-	  __cpy (&mpt1, &mpsm, p);
-	}
-      __mpsqrt (&mpsm, &mps, p);
-      mps.d[0] = X[0];
-    }
-
-  /* Evaluate a truncated power series for Atan(s).  */
-  n = __atan_np[p];
-  mptwoim1.d[1] = __atan_twonm1[p].d;
-  __dvd (&mpsm, &mptwoim1, &mpt, p);
-  for (i = n - 1; i > 1; i--)
-    {
-      mptwoim1.d[1] -= 2;
-      __dvd (&mpsm, &mptwoim1, &mpt1, p);
-      __mul (&mpsm, &mpt, &mpt2, p);
-      __sub (&mpt1, &mpt2, &mpt, p);
-    }
-  __mul (&mps, &mpt, &mpt1, p);
-  __sub (&mps, &mpt1, &mpt, p);
-
-  /* Compute Atan(x).  */
-  mptwoim1.d[1] = 1 << m;
-  __mul (&mptwoim1, &mpt, y, p);
-}
diff --git a/sysdeps/ieee754/dbl-64/mpatan.h b/sysdeps/ieee754/dbl-64/mpatan.h
deleted file mode 100644
index 5f866a77f7..0000000000
--- a/sysdeps/ieee754/dbl-64/mpatan.h
+++ /dev/null
@@ -1,145 +0,0 @@
-/*
- * IBM Accurate Mathematical Library
- * Written by International Business Machines Corp.
- * Copyright (C) 2001-2021 Free Software Foundation, Inc.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or
- * (at your option) any later version.
- *
- * This program 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program; if not, see <https://www.gnu.org/licenses/>.
- */
-
-/******************************************************************/
-/*                                                                */
-/* MODULE_NAME:mpatan.h                                           */
-/*                                                                */
-/* common data and variables prototype and definition             */
-/******************************************************************/
-
-#ifndef MPATAN_H
-#define MPATAN_H
-
-extern const number __atan_xm[8] attribute_hidden;
-extern const number __atan_twonm1[33] attribute_hidden;
-extern const number __atan_twom[8] attribute_hidden;
-extern const int __atan_np[33] attribute_hidden;
-
-
-#ifndef AVOID_MPATAN_H
-#ifdef BIG_ENDI
-  const number
-	    __atan_xm[8] = {                         /* x[m]   */
-/**/                  {{0x00000000, 0x00000000} }, /* 0.0    */
-/**/                  {{0x3f8930be, 0x00000000} }, /* 0.0123 */
-/**/                  {{0x3f991687, 0x00000000} }, /* 0.0245 */
-/**/                  {{0x3fa923a2, 0x00000000} }, /* 0.0491 */
-/**/                  {{0x3fb930be, 0x00000000} }, /* 0.0984 */
-/**/                  {{0x3fc95810, 0x00000000} }, /* 0.198  */
-/**/                  {{0x3fda7ef9, 0x00000000} }, /* 0.414  */
-/**/                  {{0x3ff00000, 0x00000000} }, /* 1.0    */
-		    };
-  const number
- __atan_twonm1[33] = {                             /* 2n-1   */
-/**/                  {{0x00000000, 0x00000000} }, /* 0      */
-/**/                  {{0x00000000, 0x00000000} }, /* 0      */
-/**/                  {{0x00000000, 0x00000000} }, /* 0      */
-/**/                  {{0x00000000, 0x00000000} }, /* 0      */
-/**/                  {{0x40260000, 0x00000000} }, /* 11     */
-/**/                  {{0x402e0000, 0x00000000} }, /* 15     */
-/**/                  {{0x40330000, 0x00000000} }, /* 19     */
-/**/                  {{0x40350000, 0x00000000} }, /* 21     */
-/**/                  {{0x40390000, 0x00000000} }, /* 25     */
-/**/                  {{0x403d0000, 0x00000000} }, /* 29     */
-/**/                  {{0x40408000, 0x00000000} }, /* 33     */
-/**/                  {{0x40428000, 0x00000000} }, /* 37     */
-/**/                  {{0x40448000, 0x00000000} }, /* 41     */
-/**/                  {{0x40468000, 0x00000000} }, /* 45     */
-/**/                  {{0x40488000, 0x00000000} }, /* 49     */
-/**/                  {{0x404a8000, 0x00000000} }, /* 53     */
-/**/                  {{0x404b8000, 0x00000000} }, /* 55     */
-/**/                  {{0x404d8000, 0x00000000} }, /* 59     */
-/**/                  {{0x404f8000, 0x00000000} }, /* 63     */
-/**/                  {{0x4050c000, 0x00000000} }, /* 67     */
-/**/                  {{0x4051c000, 0x00000000} }, /* 71     */
-/**/                  {{0x4052c000, 0x00000000} }, /* 75     */
-/**/                  {{0x4053c000, 0x00000000} }, /* 79     */
-/**/                  {{0x4054c000, 0x00000000} }, /* 83     */
-/**/                  {{0x40554000, 0x00000000} }, /* 85     */
-/**/                  {{0x40564000, 0x00000000} }, /* 89     */
-/**/                  {{0x40574000, 0x00000000} }, /* 93     */
-/**/                  {{0x40584000, 0x00000000} }, /* 97     */
-/**/                  {{0x40594000, 0x00000000} }, /* 101    */
-/**/                  {{0x405a4000, 0x00000000} }, /* 105    */
-/**/                  {{0x405b4000, 0x00000000} }, /* 109    */
-/**/                  {{0x405c4000, 0x00000000} }, /* 113    */
-/**/                  {{0x405d4000, 0x00000000} }, /* 117    */
-		    };
-
-#else
-#ifdef LITTLE_ENDI
-
-  const number
-      __atan_xm[8] = {                             /* x[m]   */
-/**/                  {{0x00000000, 0x00000000} }, /* 0.0    */
-/**/                  {{0x00000000, 0x3f8930be} }, /* 0.0123 */
-/**/                  {{0x00000000, 0x3f991687} }, /* 0.0245 */
-/**/                  {{0x00000000, 0x3fa923a2} }, /* 0.0491 */
-/**/                  {{0x00000000, 0x3fb930be} }, /* 0.0984 */
-/**/                  {{0x00000000, 0x3fc95810} }, /* 0.198  */
-/**/                  {{0x00000000, 0x3fda7ef9} }, /* 0.414  */
-/**/                  {{0x00000000, 0x3ff00000} }, /* 1.0    */
-		    };
-  const number
-__atan_twonm1[33] = {                             /* 2n-1   */
-/**/                  {{0x00000000, 0x00000000} }, /* 0      */
-/**/                  {{0x00000000, 0x00000000} }, /* 0      */
-/**/                  {{0x00000000, 0x00000000} }, /* 0      */
-/**/                  {{0x00000000, 0x00000000} }, /* 0      */
-/**/                  {{0x00000000, 0x40260000} }, /* 11     */
-/**/                  {{0x00000000, 0x402e0000} }, /* 15     */
-/**/                  {{0x00000000, 0x40330000} }, /* 19     */
-/**/                  {{0x00000000, 0x40350000} }, /* 21     */
-/**/                  {{0x00000000, 0x40390000} }, /* 25     */
-/**/                  {{0x00000000, 0x403d0000} }, /* 29     */
-/**/                  {{0x00000000, 0x40408000} }, /* 33     */
-/**/                  {{0x00000000, 0x40428000} }, /* 37     */
-/**/                  {{0x00000000, 0x40448000} }, /* 41     */
-/**/                  {{0x00000000, 0x40468000} }, /* 45     */
-/**/                  {{0x00000000, 0x40488000} }, /* 49     */
-/**/                  {{0x00000000, 0x404a8000} }, /* 53     */
-/**/                  {{0x00000000, 0x404b8000} }, /* 55     */
-/**/                  {{0x00000000, 0x404d8000} }, /* 59     */
-/**/                  {{0x00000000, 0x404f8000} }, /* 63     */
-/**/                  {{0x00000000, 0x4050c000} }, /* 67     */
-/**/                  {{0x00000000, 0x4051c000} }, /* 71     */
-/**/                  {{0x00000000, 0x4052c000} }, /* 75     */
-/**/                  {{0x00000000, 0x4053c000} }, /* 79     */
-/**/                  {{0x00000000, 0x4054c000} }, /* 83     */
-/**/                  {{0x00000000, 0x40554000} }, /* 85     */
-/**/                  {{0x00000000, 0x40564000} }, /* 89     */
-/**/                  {{0x00000000, 0x40574000} }, /* 93     */
-/**/                  {{0x00000000, 0x40584000} }, /* 97     */
-/**/                  {{0x00000000, 0x40594000} }, /* 101    */
-/**/                  {{0x00000000, 0x405a4000} }, /* 105    */
-/**/                  {{0x00000000, 0x405b4000} }, /* 109    */
-/**/                  {{0x00000000, 0x405c4000} }, /* 113    */
-/**/                  {{0x00000000, 0x405d4000} }, /* 117    */
-		    };
-
-#endif
-#endif
-
-  const int
-    __atan_np[33] = { 0, 0, 0, 0, 6, 8,10,11,13,15,17,19,21,23,25,27,28,
-		      30,32,34,36,38,40,42,43,45,47,49,51,53,55,57,59};
-
-#endif
-#endif
diff --git a/sysdeps/ieee754/dbl-64/mpatan2.c b/sysdeps/ieee754/dbl-64/mpatan2.c
deleted file mode 100644
index 68348145af..0000000000
--- a/sysdeps/ieee754/dbl-64/mpatan2.c
+++ /dev/null
@@ -1,67 +0,0 @@
-/*
- * IBM Accurate Mathematical Library
- * written by International Business Machines Corp.
- * Copyright (C) 2001-2021 Free Software Foundation, Inc.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or
- * (at your option) any later version.
- *
- * This program 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program; if not, see <https://www.gnu.org/licenses/>.
- */
-/******************************************************************/
-/*  MODULE_NAME: mpatan2.c                                        */
-/*                                                                */
-/*  FUNCTIONS:mpatan2                                             */
-/*                                                                */
-/*  FILES NEEDED: mpa.h                                           */
-/*                mpa.c mpatan.c mpsqrt.c                         */
-/*                                                                */
-/* Multi-Precision Atan2(y,x) function subroutine,                */
-/* for precision p >= 4.                                          */
-/* y=0 is not permitted if x<=0. No error messages are given.     */
-/* The relative error of the result is bounded by 44.84*r**(1-p)  */
-/* if x <= 0,  y != 0 and by 37.33*r**(1-p) if x>0. here r=2**24. */
-/*                                                                */
-/******************************************************************/
-
-#include "mpa.h"
-
-#ifndef SECTION
-# define SECTION
-#endif
-
-/* Multi-Precision Atan2 (y, x) function subroutine, for p >= 4.
-   y = 0 is not permitted if x <= 0. No error messages are given.  */
-void
-SECTION
-__mpatan2 (mp_no *y, mp_no *x, mp_no *z, int p)
-{
-  mp_no mpt1, mpt2, mpt3;
-
-  if (X[0] <= 0)
-    {
-      __dvd (x, y, &mpt1, p);
-      __mul (&mpt1, &mpt1, &mpt2, p);
-      if (mpt1.d[0] != 0)
-	mpt1.d[0] = 1;
-      __add (&mpt2, &__mpone, &mpt3, p);
-      __mpsqrt (&mpt3, &mpt2, p);
-      __add (&mpt1, &mpt2, &mpt3, p);
-      mpt3.d[0] = Y[0];
-      __mpatan (&mpt3, &mpt1, p);
-      __add (&mpt1, &mpt1, z, p);
-    }
-  else
-    {
-      __dvd (y, x, &mpt1, p);
-      __mpatan (&mpt1, z, p);
-    }
-}
diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.c b/sysdeps/ieee754/dbl-64/mpsqrt.c
deleted file mode 100644
index 2c32490c1f..0000000000
--- a/sysdeps/ieee754/dbl-64/mpsqrt.c
+++ /dev/null
@@ -1,111 +0,0 @@
-/*
- * IBM Accurate Mathematical Library
- * written by International Business Machines Corp.
- * Copyright (C) 2001-2021 Free Software Foundation, Inc.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or
- * (at your option) any later version.
- *
- * This program 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program; if not, see <https://www.gnu.org/licenses/>.
- */
-/****************************************************************************/
-/*  MODULE_NAME:mpsqrt.c                                                    */
-/*                                                                          */
-/*  FUNCTION:mpsqrt                                                         */
-/*           fastiroot                                                      */
-/*                                                                          */
-/* FILES NEEDED:endian.h mpa.h mpsqrt.h                                     */
-/*              mpa.c                                                       */
-/* Multi-Precision square root function subroutine for precision p >= 4.    */
-/* The relative error is bounded by 3.501*r**(1-p), where r=2**24.          */
-/*                                                                          */
-/****************************************************************************/
-#include "endian.h"
-#include "mpa.h"
-
-#ifndef SECTION
-# define SECTION
-#endif
-
-#include "mpsqrt.h"
-
-/****************************************************************************/
-/* Multi-Precision square root function subroutine for precision p >= 4.    */
-/* The relative error is bounded by 3.501*r**(1-p), where r=2**24.          */
-/* Routine receives two pointers to  Multi Precision numbers:               */
-/* x (left argument) and y (next argument). Routine also receives precision */
-/* p as integer. Routine computes sqrt(*x) and stores result in *y          */
-/****************************************************************************/
-
-static double fastiroot (double);
-
-void
-SECTION
-__mpsqrt (mp_no *x, mp_no *y, int p)
-{
-  int i, m, ey;
-  double dx, dy;
-  static const mp_no mphalf = {0, {1.0, HALFRAD}};
-  static const mp_no mp3halfs = {1, {1.0, 1.0, HALFRAD}};
-  mp_no mpxn, mpz, mpu, mpt1, mpt2;
-
-  ey = EX / 2;
-  __cpy (x, &mpxn, p);
-  mpxn.e -= (ey + ey);
-  __mp_dbl (&mpxn, &dx, p);
-  dy = fastiroot (dx);
-  __dbl_mp (dy, &mpu, p);
-  __mul (&mpxn, &mphalf, &mpz, p);
-
-  m = __mpsqrt_mp[p];
-  for (i = 0; i < m; i++)
-    {
-      __sqr (&mpu, &mpt1, p);
-      __mul (&mpt1, &mpz, &mpt2, p);
-      __sub (&mp3halfs, &mpt2, &mpt1, p);
-      __mul (&mpu, &mpt1, &mpt2, p);
-      __cpy (&mpt2, &mpu, p);
-    }
-  __mul (&mpxn, &mpu, y, p);
-  EY += ey;
-}
-
-/***********************************************************/
-/* Compute a double precision approximation for 1/sqrt(x)  */
-/* with the relative error bounded by 2**-51.              */
-/***********************************************************/
-static double
-SECTION
-fastiroot (double x)
-{
-  union
-  {
-    int i[2];
-    double d;
-  } p, q;
-  double y, z, t;
-  int n;
-  static const double c0 = 0.99674, c1 = -0.53380;
-  static const double c2 = 0.45472, c3 = -0.21553;
-
-  p.d = x;
-  p.i[HIGH_HALF] = (p.i[HIGH_HALF] & 0x3FFFFFFF) | 0x3FE00000;
-  q.d = x;
-  y = p.d;
-  z = y - 1.0;
-  n = (q.i[HIGH_HALF] - p.i[HIGH_HALF]) >> 1;
-  z = ((c3 * z + c2) * z + c1) * z + c0;	/* 2**-7         */
-  z = z * (1.5 - 0.5 * y * z * z);		/* 2**-14        */
-  p.d = z * (1.5 - 0.5 * y * z * z);		/* 2**-28        */
-  p.i[HIGH_HALF] -= n;
-  t = x * p.d;
-  return p.d * (1.5 - 0.5 * p.d * t);
-}
diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.h b/sysdeps/ieee754/dbl-64/mpsqrt.h
deleted file mode 100644
index d66fc99395..0000000000
--- a/sysdeps/ieee754/dbl-64/mpsqrt.h
+++ /dev/null
@@ -1,38 +0,0 @@
-/*
- * IBM Accurate Mathematical Library
- * Written by International Business Machines Corp.
- * Copyright (C) 2001-2021 Free Software Foundation, Inc.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or
- * (at your option) any later version.
- *
- * This program 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program; if not, see <https://www.gnu.org/licenses/>.
- */
-
-/******************************************************************/
-/*                                                                */
-/* MODULE_NAME:mpatan.h                                           */
-/*                                                                */
-/* common data and variables prototype and definition             */
-/******************************************************************/
-
-#ifndef MPSQRT_H
-#define MPSQRT_H
-
-extern const int __mpsqrt_mp[33] attribute_hidden;
-
-
-#ifndef AVOID_MPSQRT_H
-  const int __mpsqrt_mp[33] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,
-			     4,4,4,4,4,4,4,4,4};
-#endif
-
-#endif
diff --git a/sysdeps/ieee754/dbl-64/mptan.c b/sysdeps/ieee754/dbl-64/mptan.c
deleted file mode 100644
index 6cc3ee58a0..0000000000
--- a/sysdeps/ieee754/dbl-64/mptan.c
+++ /dev/null
@@ -1,63 +0,0 @@
-/*
- * IBM Accurate Mathematical Library
- * written by International Business Machines Corp.
- * Copyright (C) 2001-2021 Free Software Foundation, Inc.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or
- * (at your option) any later version.
- *
- * This program 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program; if not, see <https://www.gnu.org/licenses/>.
- */
-/**********************************************************************/
-/* MODULE_NAME:mptan.c                                                */
-/*                                                                    */
-/* FUNCTION: mptan                                                    */
-/*                                                                    */
-/* FILES NEEDED: endian.h  mpa.h                                      */
-/*               mpa.c  sincos32.c branred.c                          */
-/*                                                                    */
-/* Multi-Precision tan() function subroutine, for p=32.  It is based  */
-/* on the routines mpranred() and c32(). mpranred() performs range    */
-/* reduction of a double number x into a multiple precision number    */
-/* y, such that y=x-n*pi/2, abs(y)<pi/4, n=0,+-1,+-2,....  c32()      */
-/* computes both sin(y), cos(y).  tan(x) is either sin(y)/cos(y)      */
-/* or -cos(y)/sin(y). The precision of the result is of about 559     */
-/* significant bits.                                                  */
-/*                                                                    */
-/**********************************************************************/
-#include "endian.h"
-#include "mpa.h"
-
-#ifndef SECTION
-# define SECTION
-#endif
-
-void
-SECTION
-__mptan (double x, mp_no *mpy, int p)
-{
-  int n;
-  mp_no mpw, mpc, mps;
-
-  /* Negative or positive result.  */
-  n = __mpranred (x, &mpw, p) & 0x00000001;
-  /* Computing sin(x) and cos(x).  */
-  __c32 (&mpw, &mpc, &mps, p);
-  /* Second or fourth quarter of unit circle.  */
-  if (n)
-    {
-      __dvd (&mpc, &mps, mpy, p);
-      mpy->d[0] *= -1;
-    }
-  /* tan is negative in this area.  */
-  else
-    __dvd (&mps, &mpc, mpy, p);
-}
diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c
deleted file mode 100644
index 44a313ad76..0000000000
--- a/sysdeps/ieee754/dbl-64/sincos32.c
+++ /dev/null
@@ -1,307 +0,0 @@
-/*
- * IBM Accurate Mathematical Library
- * written by International Business Machines Corp.
- * Copyright (C) 2001-2021 Free Software Foundation, Inc.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or
- * (at your option) any later version.
- *
- * This program 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program; if not, see <https://www.gnu.org/licenses/>.
- */
-/****************************************************************/
-/*  MODULE_NAME: sincos32.c                                     */
-/*                                                              */
-/*  FUNCTIONS: ss32                                             */
-/*             cc32                                             */
-/*             c32                                              */
-/*             sin32                                            */
-/*             cos32                                            */
-/*             mpsin                                            */
-/*             mpcos                                            */
-/*             mpranred                                         */
-/*             mpsin1                                           */
-/*             mpcos1                                           */
-/*                                                              */
-/* FILES NEEDED: endian.h mpa.h sincos32.h                      */
-/*               mpa.c                                          */
-/*                                                              */
-/* Multi Precision sin() and cos() function with p=32  for sin()*/
-/* cos() arcsin() and arccos() routines                         */
-/* In addition mpranred() routine  performs range  reduction of */
-/* a double number x into multi precision number   y,           */
-/* such that y=x-n*pi/2, abs(y)<pi/4,  n=0,+-1,+-2,....         */
-/****************************************************************/
-#include "endian.h"
-#include "mpa.h"
-#include "sincos32.h"
-#include <math.h>
-#include <math_private.h>
-#include <stap-probe.h>
-
-#ifndef SECTION
-# define SECTION
-#endif
-
-/* Compute Multi-Precision sin() function for given p.  Receive Multi Precision
-   number x and result stored at y.  */
-static void
-SECTION
-ss32 (mp_no *x, mp_no *y, int p)
-{
-  int i;
-  double a;
-  mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}};
-  for (i = 1; i <= p; i++)
-    mpk.d[i] = 0;
-
-  __sqr (x, &x2, p);
-  __cpy (&oofac27, &gor, p);
-  __cpy (&gor, &sum, p);
-  for (a = 27.0; a > 1.0; a -= 2.0)
-    {
-      mpk.d[1] = a * (a - 1.0);
-      __mul (&gor, &mpk, &mpt1, p);
-      __cpy (&mpt1, &gor, p);
-      __mul (&x2, &sum, &mpt1, p);
-      __sub (&gor, &mpt1, &sum, p);
-    }
-  __mul (x, &sum, y, p);
-}
-
-/* Compute Multi-Precision cos() function for given p. Receive Multi Precision
-   number x and result stored at y.  */
-static void
-SECTION
-cc32 (mp_no *x, mp_no *y, int p)
-{
-  int i;
-  double a;
-  mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}};
-  for (i = 1; i <= p; i++)
-    mpk.d[i] = 0;
-
-  __sqr (x, &x2, p);
-  mpk.d[1] = 27.0;
-  __mul (&oofac27, &mpk, &gor, p);
-  __cpy (&gor, &sum, p);
-  for (a = 26.0; a > 2.0; a -= 2.0)
-    {
-      mpk.d[1] = a * (a - 1.0);
-      __mul (&gor, &mpk, &mpt1, p);
-      __cpy (&mpt1, &gor, p);
-      __mul (&x2, &sum, &mpt1, p);
-      __sub (&gor, &mpt1, &sum, p);
-    }
-  __mul (&x2, &sum, y, p);
-}
-
-/* Compute both sin(x), cos(x) as Multi precision numbers.  */
-void
-SECTION
-__c32 (mp_no *x, mp_no *y, mp_no *z, int p)
-{
-  mp_no u, t, t1, t2, c, s;
-  int i;
-  __cpy (x, &u, p);
-  u.e = u.e - 1;
-  cc32 (&u, &c, p);
-  ss32 (&u, &s, p);
-  for (i = 0; i < 24; i++)
-    {
-      __mul (&c, &s, &t, p);
-      __sub (&s, &t, &t1, p);
-      __add (&t1, &t1, &s, p);
-      __sub (&__mptwo, &c, &t1, p);
-      __mul (&t1, &c, &t2, p);
-      __add (&t2, &t2, &c, p);
-    }
-  __sub (&__mpone, &c, y, p);
-  __cpy (&s, z, p);
-}
-
-/* Compute sin() of double-length number (X + DX) as Multi Precision number and
-   return result as double.  If REDUCE_RANGE is true, X is assumed to be the
-   original input and DX is ignored.  */
-double
-SECTION
-__mpsin (double x, double dx, bool reduce_range)
-{
-  double y;
-  mp_no a, b, c, s;
-  int n;
-  int p = 32;
-
-  if (reduce_range)
-    {
-      n = __mpranred (x, &a, p);	/* n is 0, 1, 2 or 3.  */
-      __c32 (&a, &c, &s, p);
-    }
-  else
-    {
-      n = -1;
-      __dbl_mp (x, &b, p);
-      __dbl_mp (dx, &c, p);
-      __add (&b, &c, &a, p);
-      if (x > 0.8)
-        {
-          __sub (&hp, &a, &b, p);
-          __c32 (&b, &s, &c, p);
-        }
-      else
-        __c32 (&a, &c, &s, p);	/* b = sin(x+dx)  */
-    }
-
-  /* Convert result based on which quarter of unit circle y is in.  */
-  switch (n)
-    {
-    case 1:
-      __mp_dbl (&c, &y, p);
-      break;
-
-    case 3:
-      __mp_dbl (&c, &y, p);
-      y = -y;
-      break;
-
-    case 2:
-      __mp_dbl (&s, &y, p);
-      y = -y;
-      break;
-
-    /* Quadrant not set, so the result must be sin (X + DX), which is also in
-       S.  */
-    case 0:
-    default:
-      __mp_dbl (&s, &y, p);
-    }
-  LIBC_PROBE (slowsin, 3, &x, &dx, &y);
-  return y;
-}
-
-/* Compute cos() of double-length number (X + DX) as Multi Precision number and
-   return result as double.  If REDUCE_RANGE is true, X is assumed to be the
-   original input and DX is ignored.  */
-double
-SECTION
-__mpcos (double x, double dx, bool reduce_range)
-{
-  double y;
-  mp_no a, b, c, s;
-  int n;
-  int p = 32;
-
-  if (reduce_range)
-    {
-      n = __mpranred (x, &a, p);	/* n is 0, 1, 2 or 3.  */
-      __c32 (&a, &c, &s, p);
-    }
-  else
-    {
-      n = -1;
-      __dbl_mp (x, &b, p);
-      __dbl_mp (dx, &c, p);
-      __add (&b, &c, &a, p);
-      if (x > 0.8)
-        {
-          __sub (&hp, &a, &b, p);
-          __c32 (&b, &s, &c, p);
-        }
-      else
-        __c32 (&a, &c, &s, p);	/* a = cos(x+dx)     */
-    }
-
-  /* Convert result based on which quarter of unit circle y is in.  */
-  switch (n)
-    {
-    case 1:
-      __mp_dbl (&s, &y, p);
-      y = -y;
-      break;
-
-    case 3:
-      __mp_dbl (&s, &y, p);
-      break;
-
-    case 2:
-      __mp_dbl (&c, &y, p);
-      y = -y;
-      break;
-
-    /* Quadrant not set, so the result must be cos (X + DX), which is also
-       stored in C.  */
-    case 0:
-    default:
-      __mp_dbl (&c, &y, p);
-    }
-  LIBC_PROBE (slowcos, 3, &x, &dx, &y);
-  return y;
-}
-
-/* Perform range reduction of a double number x into multi precision number y,
-   such that y = x - n * pi / 2, abs (y) < pi / 4, n = 0, +-1, +-2, ...
-   Return int which indicates in which quarter of circle x is.  */
-int
-SECTION
-__mpranred (double x, mp_no *y, int p)
-{
-  number v;
-  double t, xn;
-  int i, k, n;
-  mp_no a, b, c;
-
-  if (fabs (x) < 2.8e14)
-    {
-      t = (x * hpinv.d + toint.d);
-      xn = t - toint.d;
-      v.d = t;
-      n = v.i[LOW_HALF] & 3;
-      __dbl_mp (xn, &a, p);
-      __mul (&a, &hp, &b, p);
-      __dbl_mp (x, &c, p);
-      __sub (&c, &b, y, p);
-      return n;
-    }
-  else
-    {
-      /* If x is very big more precision required.  */
-      __dbl_mp (x, &a, p);
-      a.d[0] = 1.0;
-      k = a.e - 5;
-      if (k < 0)
-	k = 0;
-      b.e = -k;
-      b.d[0] = 1.0;
-      for (i = 0; i < p; i++)
-	b.d[i + 1] = toverp[i + k];
-      __mul (&a, &b, &c, p);
-      t = c.d[c.e];
-      for (i = 1; i <= p - c.e; i++)
-	c.d[i] = c.d[i + c.e];
-      for (i = p + 1 - c.e; i <= p; i++)
-	c.d[i] = 0;
-      c.e = 0;
-      if (c.d[1] >= HALFRAD)
-	{
-	  t += 1.0;
-	  __sub (&c, &__mpone, &b, p);
-	  __mul (&b, &hp, y, p);
-	}
-      else
-	__mul (&c, &hp, y, p);
-      n = (int) t;
-      if (x < 0)
-	{
-	  y->d[0] = -y->d[0];
-	  n = -n;
-	}
-      return (n & 3);
-    }
-}
diff --git a/sysdeps/ieee754/dbl-64/sincos32.h b/sysdeps/ieee754/dbl-64/sincos32.h
deleted file mode 100644
index e112329c19..0000000000
--- a/sysdeps/ieee754/dbl-64/sincos32.h
+++ /dev/null
@@ -1,81 +0,0 @@
-/*
- * IBM Accurate Mathematical Library
- * Written by International Business Machines Corp.
- * Copyright (C) 2001-2021 Free Software Foundation, Inc.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or
- * (at your option) any later version.
- *
- * This program 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program; if not, see <https://www.gnu.org/licenses/>.
- */
-
-/******************************************************************/
-/*                                                                */
-/* MODULE_NAME:sincos32.h                                         */
-/*                                                                */
-/* common data and variables prototype and definition             */
-/******************************************************************/
-
-#ifndef SINCOS32_H
-#define SINCOS32_H
-
-#ifdef BIG_ENDI
-static const number
-/**/          hpinv = {{0x3FE45F30, 0x6DC9C883}}, /*  0.63661977236758138    */
-/**/          toint = {{0x43380000, 0x00000000}}; /*  6755399441055744       */
-
-#else
-#ifdef LITTLE_ENDI
-static const number
-/**/          hpinv = {{0x6DC9C883, 0x3FE45F30}}, /*  0.63661977236758138    */
-/**/          toint = {{0x00000000, 0x43380000}}; /*  6755399441055744       */
-
-#endif
-#endif
-
-static const mp_no
-  oofac27 = {-3,{1.0,7.0,4631664.0,12006312.0,13118056.0,6538613.0,646354.0,
-   8508025.0,9131256.0,7548776.0,2529842.0,8864927.0,660489.0,15595125.0,12777885.0,
-   11618489.0,13348664.0,5486686.0,514518.0,11275535.0,4727621.0,3575562.0,
-   13579710.0,5829745.0,7531862.0,9507898.0,6915060.0,4079264.0,1907586.0,
-   6078398.0,13789314.0,5504104.0,14136.0}},
-  pi = {1,{1.0,3.0,
-    2375530.0,8947107.0,578323.0,1673774.0,225395.0,4498441.0,3678761.0,
-    10432976.0,536314.0,10021966.0,7113029.0,2630118.0,3723283.0,7847508.0,
-    6737716.0,15273068.0,12626985.0,12044668.0,5299519.0,8705461.0,11880201.0,
-    1544726.0,14014857.0,7994139.0,13709579.0,10918111.0,11906095.0,16610011.0,
-    13638367.0,12040417.0,11529578.0,2522774.0}},
-   hp = {1,{1.0, 1.0,
-    9576373.0,4473553.0,8677769.0,9225495.0,112697.0,10637828.0,
-    10227988.0,13605096.0,268157.0,5010983.0,3556514.0,9703667.0,
-    1861641.0,12312362.0,3368858.0,7636534.0,6313492.0,14410942.0,
-    2649759.0,12741338.0,14328708.0,9160971.0,7007428.0,12385677.0,
-    15243397.0,13847663.0,14341655.0,16693613.0,15207791.0,14408816.0,
-    14153397.0,1261387.0,6110792.0,2291862.0,4181138.0,5295267.0}};
-
-static const double toverp[75] = {
-  10680707.0,  7228996.0,  1387004.0,  2578385.0, 16069853.0,
-  12639074.0,  9804092.0,  4427841.0, 16666979.0, 11263675.0,
-  12935607.0,  2387514.0,  4345298.0, 14681673.0,  3074569.0,
-  13734428.0, 16653803.0,  1880361.0, 10960616.0,  8533493.0,
-   3062596.0,  8710556.0,  7349940.0,  6258241.0,  3772886.0,
-   3769171.0,  3798172.0,  8675211.0, 12450088.0,  3874808.0,
-   9961438.0,   366607.0, 15675153.0,  9132554.0,  7151469.0,
-   3571407.0,  2607881.0, 12013382.0,  4155038.0,  6285869.0,
-   7677882.0, 13102053.0, 15825725.0,   473591.0,  9065106.0,
-  15363067.0,  6271263.0,  9264392.0,  5636912.0,  4652155.0,
-   7056368.0, 13614112.0, 10155062.0,  1944035.0,  9527646.0,
-  15080200.0,  6658437.0,  6231200.0,  6832269.0, 16767104.0,
-   5075751.0,  3212806.0,  1398474.0,  7579849.0,  6349435.0,
-  12618859.0,  4703257.0, 12806093.0, 14477321.0,  2786137.0,
-  12875403.0,  9837734.0, 14528324.0, 13719321.0,   343717.0 };
-
-#endif