Re: Exp() function reloaded

From	"io_x" <a@b.c.invalid>
Newsgroups	sci.math.num-analysis, alt.lang.asm, comp.programming
References	(1 earlier) <4F8A5BF2.17C5@mindspring.com> <4f8bc299$0$1377$4fafbaef@reader2.news.tin.it> <4f8bea40$0$1381$4fafbaef@reader1.news.tin.it> <4f8c21ad$0$1375$4fafbaef@reader1.news.tin.it> <4f8c268c$0$1378$4fafbaef@reader2.news.tin.it>
Subject	Re: Exp() function reloaded
Date	2012-04-17 09:36 +0200
Message-ID	<4f8d1c4b$0$1379$4fafbaef@reader2.news.tin.it> (permalink)
Organization	TIN.IT (http://www.tin.it)

Cross-posted to 3 groups.

Show all headers | View raw

"io_x" <a@b.c.invalid> ha scritto nel messaggio
news:4f8c268c$0$1378$4fafbaef@reader2.news.tin.it...
>
> "io_x" <a@b.c.invalid> ha scritto nel messaggio
> news:4f8c21ad$0$1375$4fafbaef@reader1.news.tin.it...
>>
>> "io_x" <a@b.c.invalid> ha scritto nel messaggio
>> news:4f8bea40$0$1381$4fafbaef@reader1.news.tin.it...
>>
>> Do you like the exp function bithread threadExp()? :)
>
> i find one couple of error in the past post.. excuse me...
> than i think is not sure the aligned of u32 type in
> "wi=(u32*)((u8*)&ip+sizeof(double)); etc"
> here is ok but there??
>
> but because it is easier for me i chose the risk in this machine
> at last...

it is not ok, the thread one is the slower one, than
if x>33 can be abs(mexp(x)-exp(x))> DBL_EPSILON and this is not good...

possible i prefer fixed point because the error in float point
seems depend from magnitude of the number in a way i not like
and not understand and not find a way to precise calculate it...


#include <stdio.h>
#include <stdint.h>

#define  u64  uint64_t
#define  u32  uint32_t
#define  u16  uint16_t
#define  u8    uint8_t
#define  i64   int64_t
#define  i32   int32_t
#define  i16   int16_t
#define  i8     int8_t

// macro for keywords
#define   R    return
#define   P    printf
#define   F    for
#define   G    goto

int    con;

long double fs_expl(long double x)
{
    long unsigned n, square;
    long double b, e, old;

    con=0;
    for (square = 0; x > 1; x /= 2) {++con;
        ++square;
    }
    while (-1 > x) {++con;
        ++square;
        x /= 2;
    }
    e = b = n = 1;
    do {++con;
        b /= n++;
        b *= x;
        e += b;
        b /= n++;
        b *= x;
        old = e;
        e += b;
    } while (e > old);
    while (square-- != 0) {++con;
        e *= e;
    }
    return e;
}

#include <math.h>
#include <float.h>

long double ke   =2.7182818284590452353602874713527;
long double maxv =LDBL_MAX/ke;
u32         maxei=0;
long double maxef=0.0;

// alcune cose [quelle giuste]
// provengono da "Satoshi Tomabechi"
// per il resto non so 100% come funziona...
// 0 per errore altrimenti ritorna il numero di cicli
u32  mexp(long double* r, long double espon)
{long double  iip, ffp,  x, t, y;
 double        ip,  fp, xx;
 u32              v, w, i;
              // una volta sola...
 if(maxei==0) // calcola esponente intero massimo+1
     {F(w=0, t=1.0; w<0xFFFFFFFF; ++w)
         {if(t>=maxv) break;
          t=t*ke;
         }
      if(w==0xFFFFFFFF)  R  0;
      maxei=w; maxef=t;
     }
 if(r==0) R   0;
 *r=0.0;
      if(espon==0.0) {*r=1.0; R   1;}
 else if(espon< 0.0) x=-espon;
 else                x= espon;
 if(x>maxei)
        {if(espon<0) {*r=0.0;      R  1;}
         else        {*r=LDBL_MAX; R  0;}
        }
 if(x>DBL_MAX)  R  0;
 xx=x; fp=modf(xx, &ip);
 F(w=0, t=1.0, v=ip, i=0; w<v; ++w)
         {++i; t=t*ke; }
 ffp=fp; iip=t;
 v=2; t=1+ffp; x=ffp;
 F(;;){++i;
       x= x*(ffp/v);
       if(x<=LDBL_EPSILON)  break;
       t+=x; ++v;
      }
 t*=iip;
 if(espon>=0) *r=    t;
 else         *r=1.0/t;
 R   i;
}

#include <windows.h>
#include <time.h>

unsigned long  __stdcall thrInteger(void*  a)
{long double  t,     ip;
 u32  w, v, i, *pu;

 t=1.0; ip=*(long double*)a; i=0;
 F(w=0, v=ip; w<v; ++w)
         {++i;
          t=t*ke;
         }
  pu=(u32*)((u8*) a + 4*sizeof(u32));
 *(long double*)a=t; *pu=i;
 R  0;
}

unsigned long  __stdcall thrFractional(void*  a)
{long double  x, t,  fp;
 u32       w, v, i, *pu;

 fp=*(long double*)a; i=0;
 v=2; t=1+fp; x=fp;
 F(;;){++i;
       x= x*(fp/v);
       if(x<=LDBL_EPSILON)  break;
       t+=x; ++v;
      }
  pu=(u32*)((u8*) a + 4*sizeof(u32));
 *(long double*)a=t; *pu=i;
 R  0;
}

u32 threadExp(long double* r, long double espon)
{long double      ip[8], fp[8], t;
 double           iip, ffp, tt;
 void             *hThr[4];
 unsigned long       IDThr;
 u32        w, i, *wi, *wf;

              // una volta sola...
 if(maxei==0) // calcola esponente intero massimo+1
     {F(w=0, t=1.0; w<0xFFFFFFFF; ++w)
         {if(t>=maxv) break;
          t=t*ke;
         }
      if(w==0xFFFFFFFF)  R  0;
      maxei=w; maxef=t;
     }
 if(r==0) R   0;
 *r=0.0;
      if(espon==0.0) {*r=1.0; R   1;}
 else if(espon< 0.0) t=-espon;
 else                t= espon;
 if(t>DBL_MAX)  R  0;
 tt=t;
 ffp=modf(tt, &iip);
 ip[0]=iip; fp[0]=ffp;
 if(iip>0xFFFFFFF||ip[0]>=maxei)
           {if(espon<0) {*r=0.0;     R  1;}
            else        {*r=LDBL_MAX; R  0;}
           }
 // un thread calcola la parte intera
 hThr[0]=CreateThread(NULL,0,thrInteger, ip,0,&IDThr);
 if(hThr[0]==0)  R  0;
 // un thread calcola la parte frazionaria
 hThr[1]=CreateThread(NULL,0,thrFractional, fp, 0,&IDThr);
 if(hThr[1]==0) {CloseHandle(hThr[0]); R  0;}
 F(i=0; i<=1; ++i)
     WaitForSingleObject(hThr[i], INFINITE);
 CloseHandle(hThr[1]); CloseHandle(hThr[0]);
 wi=(u32*)((u8*)&ip+4*sizeof(u32)); // are there aligned u32 problems?
 wf=(u32*)((u8*)&fp+4*sizeof(u32));
 t=ip[0]*fp[0];
 i=*wi+*wf;
 if(espon>=0) *r=    t;
 else         *r=1.0/t;
 R  i;
}

int  main(void)
{long double    xx, dd;
 double     d, r, x, y;
 u32        i, z, w, c;

 dd=700.1929292992; d=dd;
 x=exp(d); i=mexp(&xx,d);
 P("exp(%.20f)=%.20f, %.20Lf  i=%d\n", d, x, xx, i);
 //xx=fs_expl(dd);
 //x=xx;
 //P("exp(%.20f)=%.20f i=%d\n", d, x, con);
 srand(0);
 F(i=0, c=0; i<100000; ++i)
    {ne:
        z=rand(); w=rand();
        d=z? (double)w/z: (double) w;
        if(d>33){++c; G  ne;}
        dd=d; mexp(&xx, dd);
        y=exp(d);
        if(abs(y- (double) xx)>DBL_EPSILON)
          {P("Errore: exp(%.20f)=%.20f!=%.20Lf c=%u\n", d, y, xx, c);
           break;
          }
    }
 if(i==100000) P("all ok c=%u\n", c);
 R  0;
}

Back to comp.programming | Previous | Next — Previous in thread | Next in thread | Find similar

Thread

Exp() function reloaded hopcode <hopcode@invalid.de> - 2012-04-15 04:31 +0200
  Re: Exp() function reloaded pete <pfiland@mindspring.com> - 2012-04-15 01:26 -0400
    Re: Exp() function reloaded hopcode <hopcode@invalid.de> - 2012-04-15 12:18 +0200
      Re: Exp() function reloaded pete <pfiland@mindspring.com> - 2012-04-15 09:30 -0400
    Re: Exp() function reloaded "io_x" <a@b.c.invalid> - 2012-04-16 09:01 +0200
      Re: Exp() function reloaded "io_x" <a@b.c.invalid> - 2012-04-16 11:50 +0200
        Re: Exp() function reloaded "io_x" <a@b.c.invalid> - 2012-04-16 15:46 +0200
          Re: Exp() function reloaded "io_x" <a@b.c.invalid> - 2012-04-16 16:07 +0200
            Re: Exp() function reloaded "io_x" <a@b.c.invalid> - 2012-04-17 09:36 +0200
              Re: Exp() function reloaded phreda <pabloreda@gmail.com> - 2012-04-19 19:43 -0700
                Re: Exp() function reloaded hopcode <hopcode@invalid.de> - 2012-04-24 14:48 +0200
                Re: Exp() function reloaded phreda <pabloreda@gmail.com> - 2012-04-25 12:36 -0700
        Re: Exp() function reloaded pete <pfiland@mindspring.com> - 2012-04-16 18:12 -0400

csiph-web