POV-Ray : Newsgroups : povray.programming : Superpatch/MegaPOV function optimization Server Time
25 Dec 2024 11:17:20 EST (-0500)
  Superpatch/MegaPOV function optimization (Message 1 to 2 of 2)  
From: Lummox JR
Subject: Superpatch/MegaPOV function optimization
Date: 19 Mar 2000 14:10:38
Message: <38D52762.654A@aol.com>
Since the discussion of a raw noise function arose in povray.general,
the
topic of function optimization came up. I figure it's about time to
start
trying to hone that code past the point of theory, so I could use a lot
of
criticism.

The basic technique here is to do a two-step process: 1.) Find all
simple
expressions or function calls that can be reduced to a constant. 2.)
Find
any expression that occurs more than once, and rewrite the function to
look something like this:

precalc(expression,rest of function)

The precalc() function is internal, not open to the user. Its value is
equal to the second argument, and the first argument's value can be used
within the second by calling expr(). So precalc() acts as a sort of
cache, calculating the bigger stuff first and reusing it later.
Take this simple example:

sqrt(sqr(x)+sqr(y)+sqr(z))-sqrt(sqr(x)+sqr(y)-sqr(z))

Of that function, two expressions can be precalculated:

sqr(x)+sqr(y)
sqr(z)

The result looks like this:

precalc(sqr(z),precalc(sqr(x)+sqr(y),
  sqrt(expr(1)+expr(0))-sqrt(expr(1)-expr(0))))

It looks more complex, but it calculates faster.

The following is the code I've developed to do all that. It hasn't been
tested, and it's far from complete. Comments of any kind are welcome.

Lummox JR
------------------------
/*
 *  Function optimization:
 *
 *  Step 1: Reduce constant expressions
 *          (Still to add: Include pi, which is constant but treated
 *          differently. (Do not include clock, which can change in
another
 *          frame.) Any special constants of this sort are not included
in
 *          the number queue, so special code is needed to handle them
 *          properly. For now, just treat them as variables.)
 *  Step 2: Find any repeated expressions and pre-compute them by adding
 *          precalc() and expr() functions. (The user has no access to
these
 *          functions; they are inserted automatically.)
 */
void optimize_function(void)
{
    int tag[FUNC_MAX_STACK_LEN];
    int conststack[FUNC_MAX_STACK_LEN];
    int stacksize[FUNC_MAX_STACK_LEN];
    int i,j,k,c,i2,j2,precalc_size;
    int redoflag=0;
    int nconsts,nexpr;
    int newlen,newN;
    char *new_ops_stack;
    DBL *new_number_stack;
    DBL tmp;

    /*
     *  PHASE 1
     *
     * Optimize any constant or semi-constant expressions if possible.
     * Examples:
     *  Expression      Interpretation      Simplified
     *  (3+2)           3 2 +               5
     *  -1              1 -(unary)          -1
     *  x^2             x 2 ^               x sqr   (i.e. sqr(x))
     *  log(2)          2 log               ln(2)
     */
    /* Do not deallocate any memory to simplify; just compress */
    do {
        /* Analyze the stacks */
        for(i=k=c=0;i<func->opslen;++i) {
            j=func->ops_stack[i];
            if(j==FE_CONSTANT) conststack[i]=c++;
            else conststack[i]=-1;
            j=FuncInfo[j].Flags;
            if(FI_NUMARGS(j)>0) {
                j2=FI_NUMARGS(j);   /* numargs is args-1 for 1-3 args,
but is 4 when args=4 */
                if(j2>2) j2=3;      /* this is a 4-arg function, like
func3d. j2=args-1 */
                k-=j2;              /* Result will be args-1 lower on
stack */
                }
            else if((j&FI_NUMBER)==FI_NUMBER) ++k;  /* Toss one more
onto the stack */
            else if((j&FI_UNARY)==FI_UNARY);        /* Stack size is not
affected */
            else if((j&FI_BINARY)==FI_BINARY) --k;  /* Acts like a 2-arg
function */
            stacksize[i]=k;
            }
        /* The key loop: Find something to optimize if possible, then go
back
           up and do it all over again. If nothing is found, the routine
           ends. If something is found, i is reset and the for loop is
           broken, so it returns to the top of the do..while. */
        for(i=0;i<func->opslen-1;++i) {
            /* Keep looking until we find a constant */
            if(func->ops_stack[i]!=FE_CONSTANT) continue;
            /* Record the value of the constant, then record the
operator
               that follows */
            tmp=func->number_stack[conststack[i]];
            j=func->ops_stack[i+1];
            /* If the next operator is also a constant, keep going until
a
               non-constant is found. If it is a function, see if we
have
               enough constants in a row to calculate it. */
            if(j==FE_CONSTANT) {
                for(k=i+2;k<func->opslen-1;++k)
if(func->ops_stack[k]!=FE_CONSTANT) break;
                /* k now points to the first non-constant following i */
                j=func->ops_stack[k];
                /* How many arguments go into this function/operation?
*/
                c=stacksize[k-1]-stacksize[k]+1;
                /* If there aren't enough constants to simplify with
this function,
                 * or if this is a variable (c==0), skip ahead */
                if(c==0 || c>k-i) {i=k; continue;}
                calc_stack_top=&func->number_stack[conststack[k-1]];
                FuncInfo[j].function();
                i2=k-c; func->opslen-=c; func->numlen-=c-1;
                for(k=i2+1;k<func->opslen;++k)
func->ops_stack[k]=func->ops_stack[k+c];
                for(k=conststack[i2]+1;k<func->numlen;++k)
func->number_stack[k]=func->number_stack[k+c-1];
                i=-1; break;
                }
            /* n^2 -> sqr(n), n^3 -> cub(n), n^0.5 -> sqrt(n) */
            /* Only the exponent is a constant; the base is a variable
expression */
            if(j==FE_exp_xy) {
                if(tmp==0.5) j2=FE_sqrt;
                else if(tmp==2.0) j2=FE_sqr;
                else if(tmp==3.0) j2=FE_cub;
                else j2=-1;
                if(j2>=0) {
                    --func->numlen;
                    --func->opslen;
                    for(k=conststack[i];k<func->numlen;++k)
func->number_stack[k]=func->number_stack[k+1];
                    for(k=i+1;k<func->opslen;++k)
func->ops_stack[k]=func->ops_stack[k+1];
                    func->ops_stack[i]=j2;
                    i=-1; break;
                    }
                }
            /* Unary operators preceded by only one constant (others
caught above) */
            else if((FI_NUMARGS(FuncInfo[j]))==0 &&
(FuncInfo[j].Flags&FI_UNARY)==FI_UNARY) {
                calc_stack_top=&func->number_stack[i];
                FuncInfo[j].function();
                --func->opslen;
                for(k=i+1;k<func->opslen;++k)
func->ops_stack[k]=func->ops_stack[k+1];
                i=-1; break;
                }
            }
        if(i>=func->opslen-1) break;
        } while(1);
    /*
     *  PHASE 2
     *
     *  Pre-compute any repeated expression with 2 operators (including
     *  constants) or more. Even a 2-op expression requires calculations
     *  which benefit from precomputing, even though this will increase
op
     *  stack size by 3 and the constant stack by at least 1. Eventually
     *  pre-calculation becomes more efficient for speed than for size,
but
     *  more speed is the intended result anyway.
     */
    redoflag=1;
    /* j in this loop represents the size of the repeated expression */
    for(j=(func->opslen-1)/2,precalc_size=0;j>1;--j) {
        /* If expression has been simplified so that j>size/2, reduce j
*/
        if(j>(func->opslen-precalc_size-1)/2) continue;
        /* Search for matching expressions */
        /* i in this loop represents the first index (if any) of the
           repeated expression */
        /* i+j*2<func->opslen-precalc_size-1 because there must be room
for
           the expression to be repeated and acted on. Remember, the
last
           operator is never a constant or variable (which a complete
           expression effectively is) unless it's the only one, so a
repeated
           expression will not be the last entry before the precalcs. */
        for(i=0;(i+j*2<func->opslen-precalc_size-1);++i) {
            if(func->ops_stack[i]==FE_PRECALC) continue;
            if(redoflag) {
                /* Setup stack, constant arrays */
                for(i2=k=c=0;i<func->opslen;++i2) {
                    j=func->ops_stack[i2];
                    if(j==FE_CONSTANT) conststack[i2]=c++;
                    else conststack[i2]=-1;
                    j=FuncInfo[j].Flags;
                    if(FI_NUMARGS(j)>0) {
                        j2=FI_NUMARGS(j);
                        if(j2>2) j2=3;
                        k-=j2;
                        }
                    else if((j&FI_NUMBER)==FI_NUMBER) ++k;
                    else if((j&FI_UNARY)==FI_UNARY);
                    else if((j&FI_BINARY)==FI_BINARY) --k;
                    stacksize[i2]=k;
                    tag[i2]=0;
                    }
                redoflag=0;
                }
            /* Start only on a new expression, with something pushed to
the stack */
            if(i>0 && stacksize[i-1]>=stacksize[i]) continue;
            /* Stack must end as high as it began (a complete
expression) */
            if(stacksize[i]!=stacksize[i+j-1]) continue;
            /*
             *  Valid expression check:
             *  1.) The stack can go no lower than stacksize[i]. If it
does, then
             *      i was in the middle of another expression, not at
the
             *      beginning. Skip to this new low stack point; an
expression
             *      may begin on the next operator.
             *  2.) An expression cannot contain precalc() or expr().
Neither
             *      should be recognized as part of a repeated
expression anyway
             *      (an expr() block should already have anything that
would have
             *      repeated), but we can be safe and at the same time
use that
             *      fact to skip ahead.
             */
            for(k=i+1;k<i+j;++k) {
                /* If this is not a valid expression, skip to the point
of
                   discontinuity because a new expression may begin
there  */
                if(stacksize[k]<stacksize[i]) {i=k; break;}
                /* ignore anything with precalc() or expr() */
                /* expr() shouldn't be found in a repeated expression
anyway,
                   because the most complex expressions were already
taken
                   out, but just in case... */
                if(func->ops_stack[k]==FE_PRECALC ||
func->ops_stack[k]==FE_EXPR) {i=k; break;}
                }
            /* If not a valid expression, move on */
            if(k==i) continue;
            /* We've found a complete expression -- count constants */
            for(k=i,nconsts=0;k<i+j;++k) if(conststack[k]>=0) ++nconsts;
            /* Look for matches */
            for(k=i+j;(k+j<func->opslen);++k) {    /* k is 2nd
expression start */
                if(stacksize[k]!=stacksize[k+j-1]) continue;
                for(i2=0;i2<j;++i2) {
                    /* Compare operators */
                    if(func->ops_stack[i+i2]!=func->ops_stack[k+i2])
break;
                    /* If a constant is found, compare the constants and
be sure
                       they match exactly */
                    if(conststack[i+i2]>=0) {
                       
if(func->number_stack[conststack[i+i2]]!=func->number_stack[conststack[k+i2]])
break;
                        }
                    }
                /* If expression matches, tag it and skip to next */
                if(i2==j) {tag[i]=1; tag[k]=1; k+=j-1;}
                }
            /* If a match was found, optimize it */
            if(tag[i]) {
                newlen=func->opslen-(j-2)*nexpr+j+1;
                newN=func->numlen-(nconsts-1)*(nexpr-1)+1;
                new_ops_stack = (char *) POV_MALLOC (newlen, "function
ops");
                new_number_stack = (DBL *) POV_MALLOC (sizeof(DBL)*newN,
"function num");

                /* Increment expression constants */
                for(k=1;k<func->opslen;++k) {
                    if(func->ops_stack[k]==FE_EXPR)
func->number_stack[conststack[k-1]]+=1;
                    }
                /* Count number of times expression was found */
                for(k=nexpr=0;k<func->opslen;++k) {
                    if(tag[k]) {++nexpr; k+=j-1;}
                    }
                /* Move operators into new array */
                /* Start with new expression */
                for(k=i2=j2=0;k<j;++k) {
                    new_ops_stack[i2++]=func->ops_stack[i+k];
                    /* Adjust constant stack */
                    if(conststack[i+k]>=0)
new_number_stack[j2++]=func->number_stack[conststack[i+k]];
                    }
                /* Move in rest of function */
                for(k=0;k<func->opslen;++k) {
                    if(tag[k]) {    /* Replace expression with expr(0)
*/
                        new_ops_stack[i2++]=FE_CONSTANT;
                        new_number_stack[j2++]=0;
                        new_ops_stack[i2++]=FE_EXPR;
                        k+=j-1;
                        }
                    else {
                        new_ops_stack[i2++]=func->ops_stack[k];
                        /* Adjust constant stack */
                        if(conststack[k]>=0)
new_number_stack[j2++]=func->number_stack[conststack[k]];
                        }
                    }
                /* Add precalc() function to op stack */
                new_ops_stack[i2++]=FE_PRECALC;
                /* The precalc() function is meaningless to our
redundancy
                   check, so if we know how many there are, we can just
                   ignore the last n entries on the op stack. */
                ++precalc_size;
                /* Switch to new operator/constant lists */
                POV_FREE(func->ops_stack);
                func->ops_stack=new_ops_stack;
                func->opslen=newlen;
                POV_FREE(func->number_stack);
                func->number_stack=new_number_stack;
                func->numlen=newN;
                /* Determine new i value */
                /* Next i processed will be at i+j+2. This is because j
                   are moved to the beginning of the stack, and j are
deleted
                   and replaced by expr(n) every place the expression is
                   found. In effect, i+=j, and we then skip over
expr(n). */
                i+=j+1;
                redoflag=1;
                }
            }
        }
}


Post a reply to this message

From: Nathan Kopp
Subject: Re: Superpatch/MegaPOV function optimization
Date: 21 Mar 2000 16:26:32
Message: <38d7e908$1@news.povray.org>
Lummox JR <Lum### [at] aolcom> wrote...
> Since the discussion of a raw noise function arose in povray.general,
> the
> topic of function optimization came up. I figure it's about time to
> start
> trying to hone that code past the point of theory, so I could use a lot
> of
> criticism.
>
> The basic technique here is to do a two-step process: 1.) Find all
> simple
> expressions or function calls that can be reduced to a constant. 2.)
> Find
> any expression that occurs more than once, and rewrite the function to
> look something like this:
>
> precalc(expression,rest of function)
>
> The precalc() function is internal, not open to the user. Its value is
> equal to the second argument, and the first argument's value can be used
> within the second by calling expr(). So precalc() acts as a sort of
> cache, calculating the bigger stuff first and reusing it later.
> Take this simple example:
>
> sqrt(sqr(x)+sqr(y)+sqr(z))-sqrt(sqr(x)+sqr(y)-sqr(z))
>
> Of that function, two expressions can be precalculated:
>
> sqr(x)+sqr(y)
> sqr(z)
>
> The result looks like this:
>
> precalc(sqr(z),precalc(sqr(x)+sqr(y),
>   sqrt(expr(1)+expr(0))-sqrt(expr(1)-expr(0))))
>
> It looks more complex, but it calculates faster.
>
> The following is the code I've developed to do all that. It hasn't been
> tested, and it's far from complete. Comments of any kind are welcome.

I like your ideas.  I haven't looked at your code, but the concept seems
solid.

-Nathan


Post a reply to this message

Copyright 2003-2023 Persistence of Vision Raytracer Pty. Ltd.