Return to oh_number.rr CVS log

Up to [local] / OpenXM / src / asir-contrib / packages / src

For bigfloat.

/* $OpenXM: OpenXM/src/asir-contrib/packages/src/oh_number.rr,v 1.9 2021/09/24 04:51:02 takayama Exp $ */
module oh_number;
localf rats_prec;
localf rats;
localf to_continued_fraction;
localf continued_fraction_to_rat;
localf sign, abs;
localf deval;

def sign(X) {
    return (X==0)? 0: (X<0)? -1: 1;
}

def abs(A) {
    return A>=0? A: -A;
}

def deval(A) {
    return eval(A*exp(0));
}

def rats_prec(X,K) {
    O = cons(["rev",no],getopt());
    L = to_continued_fraction(X,K|option_list=O);
    return L[2];
}

def to_continued_fraction(X,K) {
    R=getopt(round);
    Round=0;
    if (R==floor) { 
        Round=-1;
    }else if (R==ceil) {
        Round=1;
    }
    Sign = 1;
    if (X<0) {
        Sign = -1; X = -X; 
    }
/* continued fraction (X ~= reverse(L)) */
    if (ntype(X) > 1) Bound=2^number_floor(setprec()*deval(log(10)/log(2)));
    else Bound = 2^62;
    if (X==0 || (ntype(X) && Bound*X < 1)) {
        return [Sign, [0], 0];
    }
    L = [];
    if (ntype(X)==1) {
        X = deval(X);
    }
    N1 = 1; N2 = 0;
    D1 = 0; D2 = 1; 
    Prec=10^K; // absolute error is less than 1/Prec
    for(I=0; ;I++,X=1/X) {
        Q = number_floor(X);
        L = cons(Q,L);
        X = X - Q;
        N = Q*N1 + N2; N2=N1; N1=N;
        D = Q*D1 + D2; D2=D1; D1=D;
        if ( X==0 || (ntype(X) && Bound*X < 1) || Prec < D1*D2 ) {
            break;
        }
    }
    R = (2*irem(I,2)-1)*Sign*Round;
    if (X && (!Round || R < 0) ) {
        if (length(L)>1) {
            L = cdr(L); D1=D2; N1=N2;
        }else {
            N1 = L[0] + Round; D1 = 1;
            L  = [N1];
        }
    }
    if(getopt(rev)!=no) {
        L = reverse(L);
    }
    return [Sign,L,Sign*(N1/D1)];
}

/* continued fraction to rational number */
def continued_fraction_to_rat(L) {
    Sign = 1;
    if (type(L[1])==4) { /* if L ~= [Sign, List] */
        Sign = L[0];
        L = L[1];
    }
    if(getopt(rev)!=no) {
        L = reverse(L);
    }
    for(R=car(L),L=cdr(L); L!=[]; L=cdr(L)) {
        R = car(L) + 1/R;
    }
    return Sign*R;
}

def rats(X) {
   P=getopt(prec);
   if (type(P) == 1 || P == 0) return rats_prec(X,P|option_list=getopt());
   else return rats_prec(X,14|option_list=getopt());
}

endmodule;
end$
/*
ctrl("bigfloat",1);
setprec(100);
ctrl("real_digit",32);

printf("\n")$
Z=deval(@pi);
R=rats(Z);
printf("Error=%a\n",number_abs(Z-R));

printf("\n")$
Z=0.7595;
R=rats(Z);
printf("Error=%a\n",number_abs(Z-R));

printf("\n")$
Z=deval(2^(1/2));
R=rats(Z);
printf("Error=%a\n",number_abs(Z-R));

printf("\n")$
Z=deval(3^(1/2));
R=rats(Z);
printf("Error=%a\n",number_abs(Z-R));

end;
*/