Keskustelu: Koodit: recomplex -> complexX

// *** Jouni Aro ***
// *** 28.6.2014 ***

#include <math.h>
#include <stdio.h>
#include <conio.h>
#include <complex.h>

class real2d
{
    public:

    double f[2];

    real2d(int);
    real2d(int, int);

    real2d(void);
    real2d(double);
    real2d(double *fe);
    real2d(double, double);

    friend void print(real2d);
    friend double abs(real2d);
    friend double abs(double);
    friend double sgn(real2d);
    friend real2d sqrt(real2d);

    friend real2d operator-(real2d);
    friend real2d operator+(real2d, real2d);
    friend real2d operator-(real2d, real2d);

    friend real2d operator*(real2d, real2d);
    friend real2d operator*(double, real2d);
    friend real2d operator*(real2d, double);

    friend real2d operator/(real2d, real2d);
    friend real2d operator/(double, real2d);
    friend real2d operator/(real2d, double);

    /* mgn = luvun suuruus */
    friend double mgn(real2d);
    friend real2d ln(real2d);
    friend real2d exp(real2d);
    friend real2d sin(real2d);
    friend real2d cos(real2d);
    friend real2d arcsin(real2d);
    friend real2d arccos(real2d);
    friend real2d pow(real2d, int);
    friend real2d pow(real2d, double);
    friend real2d pow(real2d, real2d);
    friend void mem0x0set(void*, int);
};

void mem0x0set(void *ptr, int n)
{
    char *c=(char*)ptr;
    for (--n; n>=0; n--)
    c[n] = (char) 0x00;
}

real2d::real2d(void)
{
    mem0x0set(this, sizeof(real2d));
}

real2d::real2d(int a)
{
    f[0]=a; f[1]=0;
}

real2d::real2d(double a)
{
    f[0]=a; f[1]=0;
}

real2d::real2d(double *fe)
{
    f[0]=fe[0]; f[1]=fe[1];
}

real2d::real2d(int a, int b)
{
    f[0]=a; f[1]=b;
}

real2d::real2d(double a, double b)
{
    f[0]=a; f[1]=b;
}

double abs(double x)
{
    return x<0.0? -x: x;
}

double abs(real2d a)
{
    return abs(a.f[0]-a.f[1]);
}

double sgn(real2d a)
{
    double x=a.f[0];
    double y=a.f[1];

    if (x>=0.0)
    {
        return abs(y)<=x? 1: -1;
    }
    else
    {
        return -1;
    }
}

double mgn(real2d a)
{
    return abs(a.f[0])+abs(a.f[1]);
}

real2d sqrt(real2d c)
{
    real2d x=c;
    for (int i=0; i<64; i++)
    x = x - (x*x-c)/(2.0*x);
    return (x);
}

real2d operator-(real2d a)
{
    for (int i=0; i<2; i++)
    a.f[i]=-a.f[i];
    return a;
}

real2d operator+(real2d a, real2d b)
{
    for (int i=0; i<2; i++)
    a.f[i]+=b.f[i];
    return a;
}

real2d operator-(real2d a, real2d b)
{
    for (int i=0; i<2; i++)
    a.f[i]-=b.f[i];
    return a;
}

real2d operator*(real2d a, real2d b)
{
    double x[4];
    double y[4];
    double t[4];

    mem0x0set(x, sizeof(double)*4);
    mem0x0set(y, sizeof(double)*4);
    mem0x0set(t, sizeof(double)*4);

    for (int i=0; i<2; i++)
    {
        x[i*2]=a.f[i];
        y[i*2]=b.f[i];
    }

    for (int i=0; i<4; i++)
    for (int j=0; j<4; j++)
    t[i^j]+=x[i]*y[j];

    for (int i=0; i<2; i++)
    t[i]=t[i*2]-t[i*2+1];
    return real2d( t );
}

real2d operator*(double k, real2d b)
{
    real2d a(k);
    return a*b;
}

real2d operator*(real2d a, double k)
{
    real2d b(k);
    return a*b;
}

real2d operator/(double k, real2d b)
{
    real2d a(k);
    return a/b;
}

real2d operator/(real2d a, double k)
{
    real2d b(k);
    return a/b;
}

real2d operator/(real2d u, real2d v)
{
    double x[4];
    double y[4];

    mem0x0set(x, sizeof(double)*4);
    mem0x0set(y, sizeof(double)*4);

    for (int i=0; i<2; i++)
    {
        x[i*2]=u.f[i];
        y[i*2]=v.f[i];
    }

    double X[4];
    double K[4*4+4];

    for (int n=0, p=0; n<4; n++)
    {
        for (int i=0; i<4; i++)
        for (int j=0; j<4; j++)
        {
            if (n == (i^j))
            {
                K[p]=y[j]; ++p;
            }
        }
        K[p]=x[n]; ++p;
    }

    // gsolve ratkaisee lineaarisen yhtälöryhmän
    // Gaussin eliminointimenetelmällä...
    int gsolve(double*, double*, int);
    gsolve(X, (double*)K, 0x4);
    for (int i=0; i<2; i++)
    X[i]=X[i*2]-X[i*2+1];
    return real2d( X );
}

////////////////////////////////////////////////
// qsolve palautusarvo:                       //
//    1 <=> yhtälöryhmälle oli ratkaisu       //
//    0 <=> yhtälöryhmälle ei ollut ratkaisua //
////////////////////////////////////////////////
int gsolve(double *X, double *f, int n)
{
    double k, J;
    double *max, *mux;
    int i, j, x, M=(int)n+1;
    void swapRows(double*, double*, int);

    for (i=0; i<n-1; i++)
    {
        max=f+i+i*M;
        for (j=i+0x01; j < n; j++)
        if (abs(*(mux=f+i+j*M))>abs(*max)) max=mux;

        if (max!=(mux=f+i+i*M))
        swapRows(max, mux, M-i);

        if (abs(J=f[i+i*M])<(double)1e-15)
        return 0; /* eipä ole ratkaisuva */

        for (j=i+1; j<n; j++)
        {
            if (f[i+j*M])
            {
                k=-f[i+j*M]/J;
                for (x=i+1; x<M; x++)
                f[x+j*M]=f[x+j*M]+k*f[x+i*M];
            }
        }
    }

    for (i=n-1; i>=0; i--)
    {
        k=(double)f[n+i*M];
        for (j=n-1; j > i; j--)
        k=k-(double)(X[j]*f[j+i*M]);

        if (abs(J=f[i+i*M])<(double)1e-15)
        return 0; /* eipä ole ratkaisua */
        else X[i]=k/J; /* onpa ratkaisu */
    }

    return 1;
}

void swapRows(double *a, double *b, int n)
{
    for (int i=0; i<n; i++)
    {
        double tmp=a[i];
        a[i] = b[i];
        b[i]=tmp;
    }
}

void print(real2d a)
{
    char sg0=a.f[0]<0.0? '-': '+';
    char sg1=a.f[1]<0.0? '-': '+';
    printf("(%c%0.10f, %c%0.10f)\n",
    sg0,abs(a.f[0]),sg1,abs(a.f[1]));
}

real2d ln(real2d c)
{
    real2d fx, dx;
    real2d z(1.0);
    real2d x(1.0);
    double r1, r2;
    r2=3141592653;

    for (int k=2; k<64; k++)
    for (int j=0; j<64; j++)
    {
        real2d xx=x;
        fx=z+x; dx=z;
        double t=1.00;

        for (int i=2; i<=k; i++)
        {
            dx=dx+i*xx/(t*=i);
            xx=xx*x; fx=fx+xx/t;
        }

        r1=r2;
        x=x-(fx-c)/dx;
        r2=mgn(x);

        if ((k<63) && (j>4) && (abs(r1-r2)<1e-15))
        {
            r2=3141592653;
            break;
        }
    }

    return x;
}

real2d exp(real2d x)
{
    double t=1.0;
    real2d sum(1.000);
    real2d xx=x; sum=sum+x;

    for (int i=2; i<64; i++)
    {
        xx = xx*x;
        sum=sum+xx/(t*=i);
    }
    return sum;
}

real2d sin(real2d x)
{
    double t=1.0;
    real2d sum(0.0);
    real2d y=x*x, xx=x;

    for (int i=1, sg=0; i<64; i+=2, sg++)
    {
        if (sg&1) sum=sum-xx/t;
        else /**/ sum=sum+xx/t;
        xx=xx*y; t*=(i+1)*(i+2);
    }
    return sum;
}

real2d cos(real2d x)
{
    double t=2.0;
    real2d sum(1.0);
    real2d y=x*x, xx=y;

    for (int i=2, sg=1; i<64; i+=2, sg++)
    {
        if (sg&1) sum=sum-xx/t;
        else /**/ sum=sum+xx/t;
        xx=xx*y; t*=(i+1)*(i+2);
    }
    return sum;
}

real2d arcsin(real2d x)
{
    real2d sum=x;
    real2d xx=x*x;

    for (int i=3; i<64; i+=2)
    {
        double ylos=1.0;
        double alas=1.0;

        for (int j=1; j<i; j++)
        {
            ylos*=j;
            alas*=(j+=1);
        }
        sum=sum+ylos*(x=x*xx)/(alas*i);
    }
    return sum;
}

real2d pow(real2d x, int n)
{
    if (n)
    {
        real2d t=x;
        int i=n<0? -n: n;
        for (--i; i; i--) t=t*x;
        return n>0x00? t: real2d(1.0)/t;
    }
    else
    {
        return real2d(1.0);
    }
}

real2d pow(real2d x, double n)
{
    return pow(x, real2d(n));
}

real2d pow(real2d a, real2d x)
{
    double t=1.00;
    real2d sum(1.000);
    real2d y=x*ln(a), xx=y;

    sum=sum+xx;

    for (int i=2; i<64; i++)
    {
        xx = xx*y;
        sum=sum+xx/(t*=i);
    }
    return sum;
}

/* ** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** ** */
/* ** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** ** */
/* ** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** ** */

class complex4d
{
    public:

    real2d r, i, j, k;

    complex4d(int);
    complex4d(void);
    complex4d(double);

    complex4d(real2d);
    complex4d(real2d, real2d);
    complex4d(real2d, real2d, real2d, real2d);

    friend void print(complex4d);
    friend double abs(complex4d);

    friend double abs(double c);
    friend complex4d sqrt(complex4d);
    friend void mem0x0set(void *ptr, int);

    friend complex4d operator-(complex4d);
    friend complex4d operator+(complex4d, complex4d);
    friend complex4d operator-(complex4d, complex4d);

    friend complex4d operator*(complex4d, complex4d);
    friend complex4d operator*(double fe, complex4d);
    friend complex4d operator*(complex4d, double fe);

    friend complex4d operator/(complex4d, complex4d);
    friend complex4d operator/(double fe, complex4d);
    friend complex4d operator/(complex4d, double fe);

    /* mgn=luvun suuruus */
    friend double mgn(complex4d c);
    friend complex4d ln(complex4d);
    friend complex4d exp(complex4d);
    friend complex4d sin(complex4d);
    friend complex4d cos(complex4d);
    friend complex4d arcsin(complex4d);
    friend complex4d arccos(complex4d);
    friend complex4d pow(complex4d, int);
    friend complex4d pow(complex4d, double);
    friend complex4d pow(complex4d, complex4d);
};

complex4d::complex4d(void)
{
    mem0x0set(this, sizeof(complex4d));
}

complex4d::complex4d(int a)
{
    r=real2d(a); i=j=k=real2d(0);
}

complex4d::complex4d(double a)
{
    r=real2d(a); i=j=k=real2d(0);
}

complex4d::complex4d(real2d a)
{
    r=a; i=j=k=real2d(0);
}

complex4d::complex4d(real2d a, real2d b)
{
    r=a; i=b; j=k=real2d(0);
}

complex4d::complex4d(real2d a, real2d b, real2d c, real2d d)
{
    r=a; i=b; j=c; k=d;
}

complex4d operator-(complex4d a)
{
    return complex4d(-a.r, -a.i, -a.j, -a.k);
}

complex4d operator+(complex4d a, complex4d b)
{
    return complex4d(a.r+b.r, a.i+b.i, a.j+b.j, a.k+b.k);
}

complex4d operator-(complex4d a, complex4d b)
{
    return complex4d(a.r-b.r, a.i-b.i, a.j-b.j, a.k-b.k);
}

complex4d operator*(double k, complex4d b)
{
    complex4d a(k);
    return a*b;
}

complex4d operator*(complex4d a, double k)
{
    complex4d b(k);
    return a*b;
}

complex4d operator/(double k, complex4d b)
{
    complex4d a(k);
    return a/b;
}

complex4d operator/(complex4d a, double k)
{
    complex4d b(k);
    return a/b;
}

complex4d operator*(complex4d a, complex4d b)
{
    complex4d c;

    c.r = real2d(+1,0)*a.r*b.r + real2d(-1,0)*a.i*b.i
        + real2d(0,+1)*a.j*b.j + real2d(0,-1)*a.k*b.k;

    c.i = real2d(+1,0)*a.r*b.i + real2d(+1,0)*a.i*b.r
        + real2d(+1,0)*a.j*b.k + real2d(+1,0)*a.k*b.j;

    c.j = real2d(+1,0)*a.r*b.j + real2d(0,-1)*a.i*b.k
        + real2d(+1,0)*a.j*b.r + real2d(0,-1)*a.k*b.i;

    c.k = real2d(+1,0)*a.r*b.k + real2d(0,+1)*a.i*b.j
        + real2d(0,+1)*a.j*b.i + real2d(+1,0)*a.k*b.r;

    return c;
}

complex4d operator/(complex4d a, complex4d b)
{
    complex4d L;
    real2d z(0);

    L=complex4d(b.r, b.i, -b.j, -b.k);
    a = a*L;
    b = b*L;

    L=complex4d(b.r, -b.i, z, z);
    a = a*L;
    b = b*L;

    a.r = a.r/b.r;
    a.i = a.i/b.r;
    a.j = a.j/b.r;
    a.k = a.k/b.r;

    return a;
}

double mgn(complex4d a)
{
    double sum=0;
    sum+=mgn(a.r);
    sum+=mgn(a.i);
    sum+=mgn(a.j);
    sum+=mgn(a.k);
    return sum;
}

double abs(complex4d a)
{
    complex4d L;
    real2d z(0);

    L=complex4d(a.r, a.i, -a.j, -a.k);
    a = a*L;

    L=complex4d(a.r, -a.i, z, z);
    a = a*L;

    a.r=pow(a.r, 0.25);

    return abs(a.r);
}

complex4d sqrt(complex4d c)
{
    complex4d x=c;
    for (int i=0; i<64; i++)
    x = x - (x*x-c)/(2.0*x);
    return x;
}

complex4d ln(complex4d c)
{
    complex4d fx, dx;
    complex4d z(1.0);
    complex4d x(1.0);
    double t, r1, r2;
    r2=3141592653589;

    for (int k=2; k<32; k++)
    for (int j=0; j<32; j++)
    {
        complex4d xx=x;
        fx=z+x; dx=z;
        t=(double)1;

        for (int i=2; i<=k; i++)
        {
            dx=dx+i*xx/(t*=i);
            xx=xx*x; fx=fx+xx/t;
        }

        r1=r2;
        x=x-(fx-c)/dx;
        r2=mgn(x);

        if ((k<31) && (j>4) && (abs(r1-r2)<1e-15))
        {
            r2=3141592653589;
            break;
        }
    }

    return x;
}

complex4d exp(complex4d x)
{
    double t=1.00;
    complex4d sum(1.00);
    complex4d xx=x; sum=sum+x;

    for (int i=2; i<32; i++)
    {
        xx = xx*x;
        sum=sum+xx/(t*=i);
    }
    return sum;
}

complex4d sin(complex4d x)
{
    double t=1.0;
    complex4d sum(0.0);
    complex4d y=x*x, xx=x;

    for (int i=1, sg=0; i<32; i+=2, sg++)
    {
        if (sg&1) sum=sum-xx/t;
        else /**/ sum=sum+xx/t;
        xx=xx*y; t*=(i+1)*(i+2);
    }
    return sum;
}

complex4d cos(complex4d x)
{
    double t=2.0;
    complex4d sum(1.0);
    complex4d y=x*x, xx=y;

    for (int i=2, sg=1; i<32; i+=2, sg++)
    {
        if (sg&1) sum=sum-xx/t;
        else /**/ sum=sum+xx/t;
        xx=xx*y; t*=(i+1)*(i+2);
    }
    return sum;
}

complex4d arcsin(complex4d x)
{
    complex4d sum=x;
    complex4d xx=x*x;

    for (int i=3; i<32; i+=2)
    {
        double ylos=1.0;
        double alas=1.0;

        for (int j=1; j<i; j++)
        {
            ylos*=j;
            alas*=(j+=1);
        }
        sum=sum+ylos*(x=x*xx)/(alas*i);
    }
    return sum;
}

complex4d pow(complex4d x, int n)
{
    if (n)
    {
        complex4d t=x;
        int i=n<0? -n: n;
        for (--i; i; i--) t=t*x;
        return n>0? t: complex4d(1.0)/t;
    }
    else
    {
        return complex4d(1.0);
    }
}

complex4d pow(complex4d x, double n)
{
    return pow(x, complex4d(n));
}

complex4d pow(complex4d a, complex4d x)
{
    double t=1.0;
    complex4d sum(1.0);
    complex4d y=x*ln(a), xx=y;

    sum=sum+xx;

    for (int i=2; i<32; i++)
    {
        xx = xx*y;
        sum=sum+xx/(t*=i);
    }
    return sum;
}

void print(complex4d a)
{
    printf("{\n");
        printf("   r = "); print(a.r);
        printf("   i = "); print(a.i);
        printf("   j = "); print(a.j);
        printf("   k = "); print(a.k);
    printf("}\n");
}

/* ** *** *** ** DEBUG ** *** *** ** DEBUG ** *** *** ** DEBUG ** *** *** ** */
/* ** *** *** ** DEBUG ** *** *** ** DEBUG ** *** *** ** DEBUG ** *** *** ** */
/* ** *** *** ** DEBUG ** *** *** ** DEBUG ** *** *** ** DEBUG ** *** *** ** */

// tuotetaan satunnaiskokonaisluku
unsigned long rnd(unsigned long max)
{
    static unsigned long R16A=1L;
    static unsigned long R16B=2L;
    R16A -= 0x012357bfL;
    R16B += 0xfedca201L;
    R16A += (R16A>>16)^(R16B<<16);
    R16B -= (R16A<<16)^(R16B>>16);
    return (R16A^R16B)%max;
}

// tuotetaan satunnaisliukuluku
double drand(void)
{
    unsigned long jak=0xfffffffb;
    double sg=rnd(0x02)? -1: 1;
    double x=(double)rnd(jak);
    return sg * x/(double)jak;
}

// tuotetaan satunnaisolio
real2d real2d_rnd(void)
{
    return real2d(drand(), drand());
}

// tuotetaan satunnaisolio
complex4d complex4d_rnd(void)
{
    real2d r(drand(), drand());
    real2d i(drand(), drand());
    real2d j(drand(), drand());
    real2d k(drand(), drand());
    return complex4d(r, i, j, k);
}

#include <stdio.h>
#include <conio.h>

#pragma argsused
int main(int kpl, char* asia[])
{
    complex4d x=complex4d_rnd();
    complex4d y=complex4d_rnd();
    complex4d z=complex4d_rnd();

    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
    printf("1) Kaikilla x, y on x+(y+z)=(x+y)+z (summan liitäntälaki)\n");
    print(x+(y+z));
    print((x+y)+z);
    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");

    printf("2) K:ssa on nolla-alkio 0 niin, että kaikilla x on x+0=x
    (summan neutraalialkio)\n");
    print(x+complex4d(0));
    print(x);
    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");

    printf("3) Kaikilla x on K:ssa vasta-alkio -x siten, että x+(-x)=0\n");
    print(x+(-x));
    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");

    printf("4) Kaikilla x, y on x + y = y + x (summan vaihdantalaki)\n");
    print(x+y);
    print(y+x);
    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");

    printf("5) Kaikilla x, y on x * (y + z) = x * y + x * z
    (osittelulaki 1)\n");
    print(x * (y + z));
    print(x * y + x * z);
    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");

    printf("6) Kaikilla x, y, z on x * (y * z) = (x * y) * z
    (tulon liitäntälaki)\n");
    print(x * (y * z));
    print((x * y) * z);
    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");

    printf("7) K:ssa on ykkösalkio 1 siten, että kaikilla x on 1 * x = x
    (tulon neutraalialkio)\n");
    print(complex4d(1)*x);
    print(x);
    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");

    printf("8) Kaikilla x paitsi 0:lla on K:ssa käänteisalkio x^(-1) siten,
    että x * x^(-1) = 1 (tulon käänteisalkio)\n");
    print(x * pow(x, -1));
    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");

    printf("9) Kaikilla x, y on x * y = y * x (tulon vaihdantalaki)\n");
    print(x * y);
    print(y * x);
    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");

    printf("sin^2 x + cos^2 x = 1\n");
    print(sin(x)*sin(x) + cos(x)*cos(x));
    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");

    printf("x=sqrt(y), z=x^2 => y=z\n");
    x = sqrt(y);
    z = x*x;
    print(y);
    print(z);
    y=complex4d_rnd();
    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");

    printf("x=ln(y), z=exp(x) => y=z\n");
    x = ln(y);
    z = exp(x);
    print(y);
    print(z);
    y=complex4d_rnd()/10.0;
    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");

    printf("x=sin(y), z=arcsin(x) => y=z\n");
    x = sin(y);
    z = arcsin(x);
    print(y);
    print(z);
    y=complex4d_rnd();
    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");

    double sg=sgn(y.k);
    if (sg==+1) printf("joo\n");
    if (sg==-1) printf("joo\n");

    getch();

    return 0;
}