// PROGRAM period
// f(x) を p 回反復計算したときの固定点を求める

#include "TrueBASIC.h"
int main();
void initial(double *r, int *p, double *epsilon, double *xleft, double *xright, 
             double *gleft, double *gright);
void bisection(double r, int p, double *xleft, double *xright, double *gleft, 
               double *gright);
double f(double x, double r, int p);

int main()
{
    int i, p;
    double epsilon, gleft, gright, r, x, xleft, xright;

    initial(&r, &p, &epsilon, &xleft, &xright, &gleft, &gright);
    do
    {
        bisection(r, p, &xleft, &xright, &gleft, &gright);
    } while(!(fabs(xleft - xright) < epsilon));
    x = 0.5*(xleft + xright);
    printf("explicit demonstration of period %d behavior\n", p);
    printf("\n");
    printf("%12d %20.13f\n", 0, x);            // 結果
    for(i = 1; i <= 2*p + 1; ++i)
    {
        x = f(x, r, 1);
        printf("%12d %20.13f\n", i, x);
    }

    return 0;
}

void initial(double *r, int *p, double *epsilon, double *xleft, double *xright, 
             double *gleft, double *gright)
{
    char done[_LBUFF_], Stmp1_[_LBUFF_];

    printf("control parameter r = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", r);
    printf("period = " );
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", p);
    printf("desired precision = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", epsilon);
    strcpy(done, "no");
    while(strcmp(done, "no") == 0)          // xleft と xright の間に 0 がくるまで実行
    {
        printf("guess for xleft = ");
        fgets(Stmp1_, _LBUFF_, stdin);
        sscanf(Stmp1_, "%lg", xleft);
        printf("guess for xright = ");
        fgets(Stmp1_, _LBUFF_, stdin);
        sscanf(Stmp1_, "%lg", xright);
        (*gleft) = f(*xleft, *r, *p) - (*xleft);
        (*gright) = f(*xright, *r, *p) - (*xright);
        if((*gleft)*(*gright) < 0)
            strcpy(done, "yes");
        else
            printf("range does not necessarily enclose a root\n");
    }
}

void bisection(double r, int p, double *xleft, double *xright, double *gleft, 
               double *gright)
{
    double gmid, xmid;

    // xleft と xright の中点
    xmid = 0.5*((*xleft) + (*xright));
    gmid = f(xmid, r, p) - xmid;
    if(gmid*(*gleft) > 0)
    {
        (*xleft) = xmid;                    // xleft の変更
        (*gleft) = gmid;
    }
    else
    {
        (*xright) = xmid;                   // xright の変更
        (*gright) = gmid;
    }
}

double f(double x, double r, int p)
// f は再帰呼び出しにより定義されている
{
    double f_, y;

    if(p > 1)
    {
        y = f(x, r, p-1);
        f_ = 4*r*y*(1-y);
    }
    else
        f_ = 4*r*x*(1-x);
    return f_;
}