[Bug libquadmath/79317] New: logq is returning double precision results

[ggoeckel at presby dot edu]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=79317

            Bug ID: 79317
           Summary: logq is returning double precision results
           Product: gcc
           Version: 5.4.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: libquadmath
          Assignee: unassigned at gcc dot gnu.org
          Reporter: ggoeckel at presby dot edu
  Target Milestone: ---

While comparing my quad-precision calculations of ln(x) with logq(x), I noticed
that the two numbers agreed only at 17 digits, leading me to assume that logq
must have a double precision bug within it.

Here is my program:

#include <sstream>
#include <string>
#include <complex>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <quadmath.h>

using namespace std;

__float128 serieslog(__float128 x)
{
 if(x==1.) return 0.; 
 x=1-x;
 __float128 sum=0.,xp=1.,t=1.;
 int n=0;
 while(fabsq(t)>1.e-33&&n<100)
 {
  n++;
  xp*=x;
  t=xp/n;
  sum+=t;
  t=fabsq(t/sum);
 }
 return -sum;
}
__float128 natlog(__float128 x)
{
 __float128
lntwo=6.9314718055994530941723212145817446e-01,fourthirds=4.,twothirds=2.;
 fourthirds/=3.;
 twothirds/=3.;
 __float128  mantissa=x;
 int exponent=0;
 while(mantissa>fourthirds)
 {
  mantissa/=1024.;
  exponent+=10;
 }
 while(mantissa<twothirds)
 {
  mantissa*=512.;
  exponent-=9;
 }
 while(mantissa>fourthirds)
 {
  mantissa/=256.;
  exponent+=8;
 }
 while(mantissa<twothirds)
 {
  mantissa*=128.;
  exponent-=7;
 }
 while(mantissa>fourthirds)
 {
  mantissa/=64.;
  exponent+=6;
 }
 while(mantissa<twothirds)
 {
  mantissa*=32.;
  exponent-=5;
 }
 while(mantissa>fourthirds)
 {
  mantissa/=16.;
  exponent+=4;
 }
 while(mantissa<twothirds)
 {
  mantissa*=8.;
  exponent-=3;
 }
 while(mantissa>fourthirds)
 {
  mantissa/=4.;
  exponent+=2;
 }
 while(mantissa<twothirds)
 {
  mantissa*=2.;
  exponent-=1;
 }
 return serieslog(mantissa)+exponent*lntwo;
}
int main()
{
 char buf[128];
 fstream
  bout;
 bout.open("logtest.txt",ios::out);
 __float128 y,z,w;
 for(__float128 x=2;x<2000002;x++)
 {
  quadmath_snprintf (buf, sizeof buf, "%#*.34Qf", 40, x);
  cout << buf << endl;
  w=logq(x);
  y=fabsq((natlog(x)-w)/w);
  if (y>0.)
  {
   y=log10q(y);
   z=floorq(y);
   y=y-z;
   y=expq(logq(10.)*y);
   if(y<5.)
    z=-z;
   else
    z=-z-1.;
  }
  else
   z=34;
  quadmath_snprintf (buf, sizeof buf, "%#*.0Qf", 2, z);
  bout << buf << endl;
 }
 return 0;
}

I am using cygwin. \cygwin64\lib\gcc\x86_64-w64-mingw32\5.4.0

My compiling commands are
g++ quadlogtest.cpp -lquadmath -o quadlogtest.exe