/* Copyright (C) 2001-2023 Artifex Software, Inc.

   All Rights Reserved.

   This software is provided AS-IS with no warranty, either express or

   implied.

   This software is distributed under license and may not be copied,

   modified or distributed except as expressly authorized under the terms

   of the license contained in the file LICENSE in this distribution.

   Refer to licensing information at http://www.artifex.com or contact

   Artifex Software, Inc.,  39 Mesa Street, Suite 108A, San Francisco,

   CA 94129, USA, for further information.

*/

/* Rendering for non-mesh shadings */

#include "math_.h"

#include "memory_.h"

#include "gx.h"

#include "gserrors.h"

#include "gsmatrix.h"   /* for gscoord.h */

#include "gscoord.h"

#include "gspath.h"

#include "gsptype2.h"

#include "gxcspace.h"

#include "gxdcolor.h"

#include "gxfarith.h"

#include "gxfixed.h"

#include "gxgstate.h"

#include "gxpath.h"

#include "gxshade.h"

#include "gxdevcli.h"

#include "gxshade4.h"

#include "gsicc_cache.h"

/* ---------------- Function-based shading ---------------- */

typedef struct Fb_frame_s { /* A rudiment of old code. */

    gs_rect region;

    gs_client_color cc[4];  /* colors at 4 corners */

    int state;

} Fb_frame_t;

typedef struct Fb_fill_state_s {

    shading_fill_state_common;

    const gs_shading_Fb_t *psh;

    gs_matrix_fixed ptm;  /* parameter space -> device space */

    Fb_frame_t frame;

} Fb_fill_state_t;

/****** NEED GC DESCRIPTOR ******/

static inline void

make_other_poles(patch_curve_t curve[4])

{

    int i, j;

    for (i = 0; i < 4; i++) {

        j = (i + 1) % 4;

        curve[i].control[0].x = (curve[i].vertex.p.x * 2 + curve[j].vertex.p.x) / 3;

        curve[i].control[0].y = (curve[i].vertex.p.y * 2 + curve[j].vertex.p.y) / 3;

        curve[i].control[1].x = (curve[i].vertex.p.x + curve[j].vertex.p.x * 2) / 3;

        curve[i].control[1].y = (curve[i].vertex.p.y + curve[j].vertex.p.y * 2) / 3;

        curve[i].straight = true;

    }

}

/* Transform a point with a fixed-point result. */

static void

gs_point_transform2fixed_clamped(const gs_matrix_fixed * pmat,

                         double x, double y, gs_fixed_point * ppt)

{

    gs_point fpt;

    gs_point_transform(x, y, (const gs_matrix *)pmat, &fpt);

    ppt->x = clamp_coord(fpt.x);

    ppt->y = clamp_coord(fpt.y);

}

static int

Fb_fill_region(Fb_fill_state_t * pfs, const gs_fixed_rect *rect)

{

    patch_fill_state_t pfs1;

    patch_curve_t curve[4];

    Fb_frame_t * fp = &pfs->frame;

    int code;

    memcpy(&pfs1, (shading_fill_state_t *)pfs, sizeof(shading_fill_state_t));

    pfs1.Function = pfs->psh->params.Function;

    code = init_patch_fill_state(&pfs1);

    if (code < 0)

        return code;

    pfs1.maybe_self_intersecting = false;

    pfs1.n_color_args = 2;

    pfs1.rect = *rect;

    gs_point_transform2fixed(&pfs->ptm, fp->region.p.x, fp->region.p.y, &curve[0].vertex.p);

    gs_point_transform2fixed(&pfs->ptm, fp->region.q.x, fp->region.p.y, &curve[1].vertex.p);

    gs_point_transform2fixed(&pfs->ptm, fp->region.q.x, fp->region.q.y, &curve[2].vertex.p);

    gs_point_transform2fixed(&pfs->ptm, fp->region.p.x, fp->region.q.y, &curve[3].vertex.p);

    make_other_poles(curve);

    curve[0].vertex.cc[0] = fp->region.p.x;   curve[0].vertex.cc[1] = fp->region.p.y;

    curve[1].vertex.cc[0] = fp->region.q.x;   curve[1].vertex.cc[1] = fp->region.p.y;

    curve[2].vertex.cc[0] = fp->region.q.x;   curve[2].vertex.cc[1] = fp->region.q.y;

    curve[3].vertex.cc[0] = fp->region.p.x;   curve[3].vertex.cc[1] = fp->region.q.y;

    code = patch_fill(&pfs1, curve, NULL, NULL);

    if (term_patch_fill_state(&pfs1))

        return_error(gs_error_unregistered); /* Must not happen. */

    return code;

}

int

gs_shading_Fb_fill_rectangle(const gs_shading_t * psh0, const gs_rect * rect,

                             const gs_fixed_rect * rect_clip,

                             gx_device * dev, gs_gstate * pgs)

{

    const gs_shading_Fb_t * const psh = (const gs_shading_Fb_t *)psh0;

    gs_matrix save_ctm;

    int xi, yi, code;

    float x[2], y[2];

    Fb_fill_state_t state;

    code = shade_init_fill_state((shading_fill_state_t *) & state, psh0, dev, pgs);

    if (code < 0)

        return code;

    state.psh = psh;

    /****** HACK FOR FIXED-POINT MATRIX MULTIPLY ******/

    gs_currentmatrix((gs_gstate *) pgs, &save_ctm);

    gs_concat((gs_gstate *) pgs, &psh->params.Matrix);

    state.ptm = pgs->ctm;

    gs_setmatrix((gs_gstate *) pgs, &save_ctm);

    /* Compute the parameter X and Y ranges. */

    {

        gs_rect pbox;

        code = gs_bbox_transform_inverse(rect, &psh->params.Matrix, &pbox);

        if (code < 0)

            return code;

        x[0] = max(pbox.p.x, psh->params.Domain[0]);

        x[1] = min(pbox.q.x, psh->params.Domain[1]);

        y[0] = max(pbox.p.y, psh->params.Domain[2]);

        y[1] = min(pbox.q.y, psh->params.Domain[3]);

    }

    if (x[0] > x[1] || y[0] > y[1]) {

        /* The region is outside the shading area. */

        if (state.icclink != NULL) gsicc_release_link(state.icclink);

        return 0;

    }

    for (xi = 0; xi < 2; ++xi)

        for (yi = 0; yi < 2; ++yi) {

            float v[2];

            v[0] = x[xi], v[1] = y[yi];

            gs_function_evaluate(psh->params.Function, v,

                                 state.frame.cc[yi * 2 + xi].paint.values);

        }

    state.frame.region.p.x = x[0];

    state.frame.region.p.y = y[0];

    state.frame.region.q.x = x[1];

    state.frame.region.q.y = y[1];

    code = Fb_fill_region(&state, rect_clip);

    if (state.icclink != NULL) gsicc_release_link(state.icclink);

    return code;

}

/* ---------------- Axial shading ---------------- */

typedef struct A_fill_state_s {

    const gs_shading_A_t *psh;

    gs_point delta;

    double length;

    double t0, t1;

    double v0, v1, u0, u1;

} A_fill_state_t;

/****** NEED GC DESCRIPTOR ******/

/* Note t0 and t1 vary over [0..1], not the Domain. */

typedef struct

{

    patch_curve_t curve[4];

    gs_point corners[4];

} corners_and_curves;

/* Ghostscript cannot possibly render any patch whose bounds aren't

 * representable in fixed's. In fact, this is a larger limit than

 * we need. We notionally have an area defined by coordinates

 * that can be represented in fixed point with at least 1 bit to

 * spare.

 *

 * Any patch that lies completely outside this region can be clipped

 * away. Any patch that isn't representable by fixed points can be

 * subdivided into 4.

 *

 * This avoids us subdividing patches huge numbers of times because

 * one side is just outside the region we will accept.

 */

#define MIN_CLIP_LIMIT ((int)(fixed2int(min_fixed)/2))

#define MAX_CLIP_LIMIT ((int)(fixed2int(max_fixed)/2))

static int not_clipped_away(const gs_point *p, const gs_fixed_rect *rect)

{

    if (p[0].x < rect->p.x &&

        p[1].x < rect->p.x &&

        p[2].x < rect->p.x &&

        p[3].x < rect->p.x)

        return 0; /* Clipped away! */

    if (p[0].x > rect->q.x &&

        p[1].x > rect->q.x &&

        p[2].x > rect->q.x &&

        p[3].x > rect->q.x)

        return 0; /* Clipped away! */

    if (p[0].y < rect->p.y &&

        p[1].y < rect->p.y &&

        p[2].y < rect->p.y &&

        p[3].y < rect->p.y)

        return 0; /* Clipped away! */

    if (p[0].y > rect->q.y &&

        p[1].y > rect->q.y &&

        p[2].y > rect->q.y &&

        p[3].y > rect->q.y)

        return 0; /* Clipped away! */

    return 1;

}

#define midpoint(a,b)\

  (arith_rshift_1(a) + arith_rshift_1(b) + (((a) | (b)) & 1))

#define quarterpoint(a,b)\

  (midpoint(a,midpoint(a,b)))

static int

subdivide_patch_fill(patch_fill_state_t *pfs, patch_curve_t c[4])

{

    fixed m0, m1;

    int v0, v1;

    int changed;

    if (pfs->rect.p.x >= pfs->rect.q.x || pfs->rect.p.y >= pfs->rect.q.y)

        return 0;

    /* On entry we have a patch:

     *   c[0].vertex  c[1].vertex

     *

     *   c[3].vertex  c[2].vertex

     *

     * Only the corners are set. The control points are not!

     *

     * BUT... in terms of spacial coords, it might be different...

     * They might be flipped on X, Y or both, giving:

     *  01 or 10 or 32 or 23

     *  32    23    01    10

     * or they might be rotated, and then flipped on X, Y or both, giving:

     *  03 or 30 or 12 or 21

     *  12    21    03    30

     */

    /* The +MIDPOINT_ACCURACY in the tests below is to allow for us finding the midpoint of [a] = z+1 and [b] = z, and getting z+1,

     * and updating [a] to be z+1, hence never actually shrinking the gap. Just accept not culling the patch as

     * much as we might. See bug 706378 for an example. */

#define MIDPOINT_ACCURACY 1

#define QUARTERPOINT_ACCURACY 3

    do {

        changed = 0;

        /* Is the whole of our patch outside the clipping rectangle? */

        /* Tempting to try to roll this into the cases below, but that

         * doesn't work because we want <= or >= here. Do X ones first. */

        if ((c[0].vertex.p.x <= pfs->rect.p.x &&

             c[1].vertex.p.x <= pfs->rect.p.x &&

             c[2].vertex.p.x <= pfs->rect.p.x &&

             c[3].vertex.p.x <= pfs->rect.p.x) ||

            (c[0].vertex.p.x >= pfs->rect.q.x &&

             c[1].vertex.p.x >= pfs->rect.q.x &&

             c[2].vertex.p.x >= pfs->rect.q.x &&

             c[3].vertex.p.x >= pfs->rect.q.x))

                return 0;

        /* First, let's try to see if we can cull the patch horizontally with the clipping

         * rectangle. */

        /* Non rotated cases first. Can we cull the left hand half? */

        if (c[0].vertex.p.x < pfs->rect.p.x && c[3].vertex.p.x < pfs->rect.p.x)

        {

            /* Check 0+3 off left. */

            v0 = 0;

            v1 = 3;

            goto check_left;

        }

        else if (c[1].vertex.p.x < pfs->rect.p.x && c[2].vertex.p.x < pfs->rect.p.x)

        {

            /* Check 1+2 off left. */

            v0 = 1;

            v1 = 2;

check_left:

            /* At this point we know that the condition for the following loop is true, so it

             * can be a do...while rather than a while. */

            do

            {

                /* Let's form (X coords only):

                 *

                 * c[v0].vertex  m0  c[v0^1].vertex

                 * c[v1].vertex  m1  c[v1^1].vertex

                 */

                m0 = midpoint(c[0].vertex.p.x, c[1].vertex.p.x);

                if (m0 >= pfs->rect.p.x)

                    goto check_left_quarter;

                m1 = midpoint(c[3].vertex.p.x, c[2].vertex.p.x);

                if (m1 >= pfs->rect.p.x)

                    goto check_left_quarter;

                /* So, we can completely discard the left hand half of the patch. */

                c[v0].vertex.p.x = m0;

                c[v0].vertex.p.y = midpoint(c[0].vertex.p.y, c[1].vertex.p.y);

                c[v1].vertex.p.x = m1;

                c[v1].vertex.p.y = midpoint(c[3].vertex.p.y, c[2].vertex.p.y);

                c[v0].vertex.cc[0] = (c[0].vertex.cc[0] + c[1].vertex.cc[0])/2;

                c[v1].vertex.cc[0] = (c[3].vertex.cc[0] + c[2].vertex.cc[0])/2;

                changed = 1;

            }

            while (c[v0].vertex.p.x < pfs->rect.p.x && c[v1].vertex.p.x < pfs->rect.p.x);

            if (0)

            {

check_left_quarter:

                /* At this point we know that the condition for the following loop is true, so it

                 * can be a do...while rather than a while. */

                do

                {

                    /* Let's form (X coords only):

                     *

                     * c[v0].vertex  m0  x  x  c[v0^1].vertex

                     * c[v1].vertex  m1  x  x  c[v1^1].vertex

                     */

                    m0 = quarterpoint(c[v0].vertex.p.x, c[v0^1].vertex.p.x);

                    if (m0 >= pfs->rect.p.x)

                        break;

                    m1 = quarterpoint(c[v1].vertex.p.x, c[v1^1].vertex.p.x);

                    if (m1 >= pfs->rect.p.x)

                        break;

                    /* So, we can completely discard the left hand quarter of the patch. */

                    c[v0].vertex.p.x = m0;

                    c[v0].vertex.p.y = midpoint(c[v0].vertex.p.y, c[v0^1].vertex.p.y);

                    c[v1].vertex.p.x = m1;

                    c[v1].vertex.p.y = midpoint(c[v1].vertex.p.y, c[v1^1].vertex.p.y);

                    c[v0].vertex.cc[0] = (c[v0].vertex.cc[0] + 3*c[v0^1].vertex.cc[0])/4;

                    c[v1].vertex.cc[0] = (c[v1].vertex.cc[0] + 3*c[v1^1].vertex.cc[0])/4;

                    changed = 1;

                }

                while (c[v0].vertex.p.x < pfs->rect.p.x && c[v1].vertex.p.x < pfs->rect.p.x);

            }

        }

        /* or the right hand half? */

        if (c[0].vertex.p.x > pfs->rect.q.x && c[3].vertex.p.x > pfs->rect.q.x)

        {

            /* Check 0+3 off right. */

            v0 = 0;

            v1 = 3;

            goto check_right;

        }

        else if (c[1].vertex.p.x > pfs->rect.q.x && c[2].vertex.p.x > pfs->rect.q.x)

        {

            /* Check 1+2 off right. */

            v0 = 1;

            v1 = 2;

check_right:

            /* At this point we know that the condition for the following loop is true, so it

             * can be a do...while rather than a while. */

            do

            {

                /* Let's form (X coords only):

                 *

                 * c[v0].vertex  m0  c[v0^1].vertex

                 * c[v1].vertex  m1  c[v1^1].vertex

                 */

                m0 = midpoint(c[0].vertex.p.x, c[1].vertex.p.x);

                if (m0 <= pfs->rect.q.x+MIDPOINT_ACCURACY)

                    goto check_right_quarter;

                m1 = midpoint(c[3].vertex.p.x, c[2].vertex.p.x);

                if (m1 <= pfs->rect.q.x+MIDPOINT_ACCURACY)

                    goto check_right_quarter;

                /* So, we can completely discard the left hand half of the patch. */

                c[v0].vertex.p.x = m0;

                c[v0].vertex.p.y = midpoint(c[0].vertex.p.y, c[1].vertex.p.y);

                c[v1].vertex.p.x = m1;

                c[v1].vertex.p.y = midpoint(c[3].vertex.p.y, c[2].vertex.p.y);

                c[v0].vertex.cc[0] = (c[0].vertex.cc[0] + c[1].vertex.cc[0])/2;

                c[v1].vertex.cc[0] = (c[3].vertex.cc[0] + c[2].vertex.cc[0])/2;

                changed = 1;

            }

            while (c[v0].vertex.p.x > pfs->rect.q.x+MIDPOINT_ACCURACY && c[v1].vertex.p.x > pfs->rect.q.x+MIDPOINT_ACCURACY);

            if (0)

            {

check_right_quarter:

                /* At this point we know that the condition for the following loop is true, so it

                 * can be a do...while rather than a while. */

                do

                {

                    /* Let's form (X coords only):

                     *

                     * c[v0].vertex  m0  x  x  c[v0^1].vertex

                     * c[v1].vertex  m1  x  x  c[v1^1].vertex

                     */

                    m0 = quarterpoint(c[v0].vertex.p.x, c[v0^1].vertex.p.x);

                    if (m0 <= pfs->rect.q.x+QUARTERPOINT_ACCURACY)

                        break;

                    m1 = quarterpoint(c[v1].vertex.p.x, c[v1^1].vertex.p.x);

                    if (m1 <= pfs->rect.q.x+QUARTERPOINT_ACCURACY)

                        break;

                    /* So, we can completely discard the left hand half of the patch. */

                    c[v0].vertex.p.x = m0;

                    c[v0].vertex.p.y = quarterpoint(c[v0].vertex.p.y, c[v0^1].vertex.p.y);

                    c[v1].vertex.p.x = m1;

                    c[v1].vertex.p.y = quarterpoint(c[v1].vertex.p.y, c[v1^1].vertex.p.y);

                    c[v0].vertex.cc[0] = (c[v0].vertex.cc[0] + 3*c[v0^1].vertex.cc[0])/4;

                    c[v1].vertex.cc[0] = (c[v1].vertex.cc[0] + 3*c[v1^1].vertex.cc[0])/4;

                    changed = 1;

                }

                while (c[v0].vertex.p.x > pfs->rect.q.x+QUARTERPOINT_ACCURACY && c[v1].vertex.p.x > pfs->rect.q.x+QUARTERPOINT_ACCURACY);

            }

        }

        /* Now, rotated cases: Can we cull the left hand half? */

        if (c[0].vertex.p.x < pfs->rect.p.x && c[1].vertex.p.x < pfs->rect.p.x)

        {

            /* Check 0+1 off left. */

            v0 = 0;

            v1 = 1;

            goto check_rot_left;

        }

        else if (c[3].vertex.p.x < pfs->rect.p.x && c[2].vertex.p.x < pfs->rect.p.x)

        {

            /* Check 3+2 off left. */

            v0 = 3;

            v1 = 2;

check_rot_left:

            /* At this point we know that the condition for the following loop is true, so it

             * can be a do...while rather than a while. */

            do

            {

                /* Let's form (X coords only):

                 *

                 * c[v0].vertex    m0  c[v0^3].vertex

                 * c[v1^3].vertex  m1  c[v1].vertex

                 */

                m0 = midpoint(c[0].vertex.p.x, c[3].vertex.p.x);

                if (m0 >= pfs->rect.p.x)

                    goto check_rot_left_quarter;

                m1 = midpoint(c[1].vertex.p.x, c[2].vertex.p.x);

                if (m1 >= pfs->rect.p.x)

                    goto check_rot_left_quarter;

                /* So, we can completely discard the left hand half of the patch. */

                c[v0].vertex.p.x = m0;

                c[v0].vertex.p.y = midpoint(c[0].vertex.p.y, c[3].vertex.p.y);

                c[v1].vertex.p.x = m1;

                c[v1].vertex.p.y = midpoint(c[1].vertex.p.y, c[2].vertex.p.y);

                c[v0].vertex.cc[0] = (c[0].vertex.cc[0] + c[3].vertex.cc[0])/2;

                c[v1].vertex.cc[0] = (c[1].vertex.cc[0] + c[2].vertex.cc[0])/2;

                changed = 1;

            }

            while (c[v0].vertex.p.x < pfs->rect.p.x && c[v1].vertex.p.x < pfs->rect.p.x);

            if (0)

            {

check_rot_left_quarter:

                /* At this point we know that the condition for the following loop is true, so it

                 * can be a do...while rather than a while. */

                do

                {

                    /* Let's form (X coords only):

                     *

                     * c[v0].vertex  m0  x  x  c[v0^3].vertex

                     * c[v1].vertex  m1  x  x  c[v1^3].vertex

                     */

                    m0 = quarterpoint(c[v0].vertex.p.x, c[v0^3].vertex.p.x);

                    if (m0 >= pfs->rect.p.x)

                        break;

                    m1 = quarterpoint(c[v1].vertex.p.x, c[v1^3].vertex.p.x);

                    if (m1 >= pfs->rect.p.x)

                        break;

                    /* So, we can completely discard the left hand half of the patch. */

                    c[v0].vertex.p.x = m0;

                    c[v0].vertex.p.y = quarterpoint(c[v0].vertex.p.y, c[v0^3].vertex.p.y);

                    c[v1].vertex.p.x = m1;

                    c[v1].vertex.p.y = quarterpoint(c[v1].vertex.p.y, c[v1^3].vertex.p.y);

                    c[v0].vertex.cc[0] = (c[v0].vertex.cc[0] + 3*c[v0^3].vertex.cc[0])/4;

                    c[v1].vertex.cc[0] = (c[v1].vertex.cc[0] + 3*c[v1^3].vertex.cc[0])/4;

                    changed = 1;

                }

                while (c[v0].vertex.p.x < pfs->rect.p.x && c[v1].vertex.p.x < pfs->rect.p.x);

            }

        }

        /* or the right hand half? */

        if (c[0].vertex.p.x > pfs->rect.q.x && c[1].vertex.p.x > pfs->rect.q.x)

        {

            /* Check 0+1 off right. */

            v0 = 0;

            v1 = 1;

            goto check_rot_right;

        }

        else if (c[3].vertex.p.x > pfs->rect.q.x && c[2].vertex.p.x > pfs->rect.q.x)

        {

            /* Check 3+2 off right. */

            v0 = 3;

            v1 = 2;

check_rot_right:

            /* At this point we know that the condition for the following loop is true, so it

             * can be a do...while rather than a while. */

            do

            {

                /* Let's form (X coords only):

                 *

                 * c[v0].vertex  m0  c[v0^3].vertex

                 * c[v1].vertex  m1  c[v1^3].vertex

                 */

                m0 = midpoint(c[0].vertex.p.x, c[3].vertex.p.x);

                if (m0 <= pfs->rect.q.x+MIDPOINT_ACCURACY)

                    goto check_rot_right_quarter;

                m1 = midpoint(c[1].vertex.p.x, c[2].vertex.p.x);

                if (m1 <= pfs->rect.q.x+MIDPOINT_ACCURACY)

                    goto check_rot_right_quarter;

                /* So, we can completely discard the left hand half of the patch. */

                c[v0].vertex.p.x = m0;

                c[v0].vertex.p.y = midpoint(c[0].vertex.p.y, c[3].vertex.p.y);

                c[v1].vertex.p.x = m1;

                c[v1].vertex.p.y = midpoint(c[1].vertex.p.y, c[2].vertex.p.y);

                c[v0].vertex.cc[0] = (c[0].vertex.cc[0] + c[3].vertex.cc[0])/2;

                c[v1].vertex.cc[0] = (c[1].vertex.cc[0] + c[2].vertex.cc[0])/2;

                changed = 1;

            }

            while (c[v0].vertex.p.x > pfs->rect.q.x+MIDPOINT_ACCURACY && c[v1].vertex.p.x > pfs->rect.q.x+MIDPOINT_ACCURACY);

            if (0)

            {

check_rot_right_quarter:

                /* At this point we know that the condition for the following loop is true, so it

                 * can be a do...while rather than a while. */

                do

                {

                    /* Let's form (X coords only):

                     *

                     * c[v0].vertex  m0  c[v0^3].vertex

                     * c[v1].vertex  m1  c[v1^3].vertex

                     */

                    m0 = quarterpoint(c[v0].vertex.p.x, c[v0^3].vertex.p.x);

                    if (m0 <= pfs->rect.q.x+QUARTERPOINT_ACCURACY)

                        break;

                    m1 = quarterpoint(c[v1].vertex.p.x, c[v1^3].vertex.p.x);

                    if (m1 <= pfs->rect.q.x+QUARTERPOINT_ACCURACY)

                        break;

                    /* So, we can completely discard the left hand half of the patch. */

                    c[v0].vertex.p.x = m0;

                    c[v0].vertex.p.y = quarterpoint(c[v0].vertex.p.y, c[v0^3].vertex.p.y);

                    c[v1].vertex.p.x = m1;

                    c[v1].vertex.p.y = quarterpoint(c[v1].vertex.p.y, c[v1^3].vertex.p.y);

                    c[v0].vertex.cc[0] = (c[v0].vertex.cc[0] + 3*c[v0^3].vertex.cc[0])/4;

                    c[v1].vertex.cc[0] = (c[v1].vertex.cc[0] + 3*c[v1^3].vertex.cc[0])/4;

                    changed = 1;

                }

                while (c[v0].vertex.p.x > pfs->rect.q.x+QUARTERPOINT_ACCURACY && c[v1].vertex.p.x > pfs->rect.q.x+QUARTERPOINT_ACCURACY);

            }

        }

        /* Is the whole of our patch outside the clipping rectangle? */

        /* Tempting to try to roll this into the cases below, but that

         * doesn't work because we want <= or >= here. Do Y ones. Can't have

         * done this earlier, as the previous set of tests might have reduced

         * the range here. */

        if ((c[0].vertex.p.y <= pfs->rect.p.y &&

             c[1].vertex.p.y <= pfs->rect.p.y &&

             c[2].vertex.p.y <= pfs->rect.p.y &&

             c[3].vertex.p.y <= pfs->rect.p.y) ||

            (c[0].vertex.p.y >= pfs->rect.q.y &&

             c[1].vertex.p.y >= pfs->rect.q.y &&

             c[2].vertex.p.y >= pfs->rect.q.y &&

             c[3].vertex.p.y >= pfs->rect.q.y))

            return 0;

        /* Now, let's try to see if we can cull the patch vertically with the clipping

         * rectangle. */

        /* Non rotated cases first. Can we cull the top half? */

        if (c[0].vertex.p.y < pfs->rect.p.y && c[1].vertex.p.y < pfs->rect.p.y)

        {

            /* Check 0+1 off above. */

            v0 = 0;

            v1 = 1;

            goto check_above;

        }

        else if (c[3].vertex.p.y < pfs->rect.p.y && c[2].vertex.p.y < pfs->rect.p.y)

        {

            /* Check 3+2 off above. */

            v0 = 3;

            v1 = 2;

check_above:

            /* At this point we know that the condition for the following loop is true, so it

             * can be a do...while rather than a while. */

            do

            {

                /* Let's form (Y coords only):

                 *

                 * c[v0].vertex     c[v1].vertex

                 * m0               m1

                 * c[v0^3].vertex   c[v1^3].vertex

                 */

                m0 = midpoint(c[0].vertex.p.y, c[3].vertex.p.y);

                if (m0 >= pfs->rect.p.y)

                    goto check_above_quarter;

                m1 = midpoint(c[1].vertex.p.y, c[2].vertex.p.y);

                if (m1 >= pfs->rect.p.y)

                    goto check_above_quarter;

                /* So, we can completely discard the top half of the patch. */

                c[v0].vertex.p.x = midpoint(c[0].vertex.p.x, c[3].vertex.p.x);

                c[v0].vertex.p.y = m0;

                c[v1].vertex.p.x = midpoint(c[1].vertex.p.x, c[2].vertex.p.x);

                c[v1].vertex.p.y = m1;

                c[v0].vertex.cc[0] = (c[0].vertex.cc[0] + c[3].vertex.cc[0])/2;

                c[v1].vertex.cc[0] = (c[1].vertex.cc[0] + c[2].vertex.cc[0])/2;

                changed = 1;

            }

            while (c[v0].vertex.p.y < pfs->rect.p.y && c[v1].vertex.p.y < pfs->rect.p.y);

            if (0)

            {

check_above_quarter:

                /* At this point we know that the condition for the following loop is true, so it

                 * can be a do...while rather than a while. */

                do

                {

                    /* Let's form (Y coords only):

                     *

                     * c[v0].vertex     c[v1].vertex

                     * m0               m1

                     * x                x

                     * x                x

                     * c[v0^3].vertex   c[v1^3].vertex

                     */

                    m0 = quarterpoint(c[v0].vertex.p.y, c[v0^3].vertex.p.y);

                    if (m0 >= pfs->rect.p.y)

                        break;

                    m1 = quarterpoint(c[v1].vertex.p.y, c[v1^3].vertex.p.y);

                    if (m1 >= pfs->rect.p.y)

                        break;

                    /* So, we can completely discard the top half of the patch. */

                    c[v0].vertex.p.x = quarterpoint(c[v0].vertex.p.x, c[v0^3].vertex.p.x);

                    c[v0].vertex.p.y = m0;

                    c[v1].vertex.p.x = quarterpoint(c[v1].vertex.p.x, c[v1^3].vertex.p.x);

                    c[v1].vertex.p.y = m1;

                    c[v0].vertex.cc[0] = (c[v0].vertex.cc[0] + 3*c[v0^3].vertex.cc[0])/4;

                    c[v1].vertex.cc[0] = (c[v1].vertex.cc[0] + 3*c[v1^3].vertex.cc[0])/4;

                    changed = 1;

                }

                while (c[v0].vertex.p.y < pfs->rect.p.y && c[v1].vertex.p.y < pfs->rect.p.y);

            }

        }

        /* or the bottom half? */

        if (c[0].vertex.p.y > pfs->rect.q.y && c[1].vertex.p.y > pfs->rect.q.y)

        {

            /* Check 0+1 off bottom. */

            v0 = 0;

            v1 = 1;

            goto check_bottom;

        }

        else if (c[3].vertex.p.y > pfs->rect.q.y && c[2].vertex.p.y > pfs->rect.q.y)

        {

            /* Check 3+2 off bottom. */

            v0 = 3;

            v1 = 2;

check_bottom:

            /* At this point we know that the condition for the following loop is true, so it

             * can be a do...while rather than a while. */

            do

            {

                /* Let's form (Y coords only):

                 *

                 * c[v0].vertex     c[v1].vertex

                 * m0               m1

                 * c[v0^3].vertex   c[v1^3].vertex

                 */

                m0 = midpoint(c[0].vertex.p.y, c[3].vertex.p.y);

                if (m0 <= pfs->rect.q.y+MIDPOINT_ACCURACY)

                    goto check_bottom_quarter;

                m1 = midpoint(c[1].vertex.p.y, c[2].vertex.p.y);

                if (m1 <= pfs->rect.q.y+MIDPOINT_ACCURACY)

                    goto check_bottom_quarter;

                /* So, we can completely discard the bottom half of the patch. */

                c[v0].vertex.p.x = midpoint(c[0].vertex.p.x, c[3].vertex.p.x);

                c[v0].vertex.p.y = m0;

                c[v1].vertex.p.x = midpoint(c[1].vertex.p.x, c[2].vertex.p.x);

                c[v1].vertex.p.y = m1;

                c[v0].vertex.cc[0] = (c[0].vertex.cc[0] + c[3].vertex.cc[0])/2;

                c[v1].vertex.cc[0] = (c[1].vertex.cc[0] + c[2].vertex.cc[0])/2;

                changed = 1;

            }

            while (c[v0].vertex.p.y > pfs->rect.q.y+MIDPOINT_ACCURACY && c[v1].vertex.p.y > pfs->rect.q.y+MIDPOINT_ACCURACY);

            if (0)

            {

check_bottom_quarter:

                /* At this point we know that the condition for the following loop is true, so it

                 * can be a do...while rather than a while. */

                do

                {

                    /* Let's form (Y coords only):

                     *

                     * c[v0].vertex     c[v1].vertex

                     * x                x

                     * x                x

                     * m0               m1

                     * c[v0^3].vertex   c[v1^3].vertex

                     */

                    m0 = quarterpoint(c[v0].vertex.p.y, c[v0^3].vertex.p.y);

                    if (m0 <= pfs->rect.q.y+QUARTERPOINT_ACCURACY)

                        break;

                    m1 = quarterpoint(c[v1].vertex.p.y, c[v1^3].vertex.p.y);

                    if (m1 <= pfs->rect.q.y+QUARTERPOINT_ACCURACY)

                        break;

                    /* So, we can completely discard the bottom half of the patch. */

                    c[v0].vertex.p.x = quarterpoint(c[v0].vertex.p.x, c[v0^3].vertex.p.x);

                    c[v0].vertex.p.y = m0;

                    c[v1].vertex.p.x = quarterpoint(c[v1].vertex.p.x, c[v1^3].vertex.p.x);

                    c[v1].vertex.p.y = m1;

                    c[v0].vertex.cc[0] = (c[v0].vertex.cc[0] + 3*c[v0^3].vertex.cc[0])/4;

                    c[v1].vertex.cc[0] = (c[v1].vertex.cc[0] + 3*c[v1^3].vertex.cc[0])/4;

                    changed = 1;

                }

                while (c[v0].vertex.p.y > pfs->rect.q.y+QUARTERPOINT_ACCURACY && c[v1].vertex.p.y > pfs->rect.q.y+QUARTERPOINT_ACCURACY);

            }

        }

        /* Now, rotated cases: Can we cull the top half? */

        if (c[0].vertex.p.y < pfs->rect.p.y && c[3].vertex.p.y < pfs->rect.p.y)

        {

            /* Check 0+3 off above. */

            v0 = 0;

            v1 = 3;

            goto check_rot_above;

        }

        else if (c[1].vertex.p.y < pfs->rect.p.y && c[2].vertex.p.y < pfs->rect.p.y)

        {

            /* Check 1+2 off above. */

            v0 = 1;

            v1 = 2;

check_rot_above:

            /* At this point we know that the condition for the following loop is true, so it

             * can be a do...while rather than a while. */

            do

            {

                /* Let's form (Y coords only):

                 *

                 * c[v0].vertex     c[v1].vertex

                 * m0               m1

                 * c[v0^1].vertex   c[v1^1].vertex

                 */

                m0 = midpoint(c[0].vertex.p.y, c[1].vertex.p.y);

                if (m0 >= pfs->rect.p.y)

                    goto check_rot_above_quarter;

                m1 = midpoint(c[3].vertex.p.y, c[2].vertex.p.y);

                if (m1 >= pfs->rect.p.y)

                    goto check_rot_above_quarter;

                /* So, we can completely discard the top half of the patch. */

                c[v0].vertex.p.x = midpoint(c[0].vertex.p.x, c[1].vertex.p.x);

                c[v0].vertex.p.y = m0;

                c[v1].vertex.p.x = midpoint(c[3].vertex.p.x, c[2].vertex.p.x);

                c[v1].vertex.p.y = m1;

                c[v0].vertex.cc[0] = (c[0].vertex.cc[0] + c[1].vertex.cc[0])/2;

                c[v1].vertex.cc[0] = (c[3].vertex.cc[0] + c[2].vertex.cc[0])/2;

                changed = 1;

            }

            while (c[v0].vertex.p.y < pfs->rect.p.y && c[v1].vertex.p.y < pfs->rect.p.y);

            if (0)

            {

check_rot_above_quarter:

                /* At this point we know that the condition for the following loop is true, so it

                 * can be a do...while rather than a while. */

                do

                {

                    /* Let's form (Y coords only):

                     *

                     * c[v0].vertex     c[v1].vertex

                     * m0               m1

                     * x                x

                     * x                x

                     * c[v0^1].vertex   c[v1^1].vertex

                     */

                    m0 = quarterpoint(c[v0].vertex.p.y, c[v0^1].vertex.p.y);

                    if (m0 >= pfs->rect.p.y)

                        break;

                    m1 = quarterpoint(c[v1].vertex.p.y, c[v1^1].vertex.p.y);

                    if (m1 >= pfs->rect.p.y)

                        break;

                    /* So, we can completely discard the top half of the patch. */

                    c[v0].vertex.p.x = quarterpoint(c[v0].vertex.p.x, c[v0^1].vertex.p.x);

                    c[v0].vertex.p.y = m0;

                    c[v1].vertex.p.x = quarterpoint(c[v1].vertex.p.x, c[v1^1].vertex.p.x);

                    c[v1].vertex.p.y = m1;

                    c[v0].vertex.cc[0] = (c[v0].vertex.cc[0] + 3*c[v0^1].vertex.cc[0])/4;

                    c[v1].vertex.cc[0] = (c[v1].vertex.cc[0] + 3*c[v1^1].vertex.cc[0])/4;

                    changed = 1;

                }

                while (c[v0].vertex.p.y < pfs->rect.p.y && c[v1].vertex.p.y < pfs->rect.p.y);

            }

        }

        /* or the bottom half? */

        if (c[0].vertex.p.y > pfs->rect.q.y && c[3].vertex.p.y > pfs->rect.q.y)

        {

            /* Check 0+3 off the bottom. */

            v0 = 0;

            v1 = 3;

            goto check_rot_bottom;

        }

        else if (c[1].vertex.p.y > pfs->rect.q.y && c[2].vertex.p.y > pfs->rect.q.y)

        {

            /* Check 1+2 off the bottom. */

            v0 = 1;

            v1 = 2;

check_rot_bottom:

            /* At this point we know that the condition for the following loop is true, so it

             * can be a do...while rather than a while. */

            do

            {

                /* Let's form (Y coords only):

                 *

                 * c[v0].vertex     c[v1].vertex

                 * m0               m1

                 * c[v0^1].vertex   c[v1^1].vertex

                 */

                m0 = midpoint(c[0].vertex.p.y, c[1].vertex.p.y);

                if (m0 <= pfs->rect.q.y+MIDPOINT_ACCURACY)

                    goto check_rot_bottom_quarter;

                m1 = midpoint(c[3].vertex.p.y, c[2].vertex.p.y);

                if (m1 <= pfs->rect.q.y+MIDPOINT_ACCURACY)

                    goto check_rot_bottom_quarter;

                /* So, we can completely discard the left hand half of the patch. */

                c[v0].vertex.p.x = midpoint(c[0].vertex.p.x, c[1].vertex.p.x);

                c[v0].vertex.p.y = m0;

                c[v1].vertex.p.x = midpoint(c[3].vertex.p.x, c[2].vertex.p.x);

                c[v1].vertex.p.y = m1;

                c[v0].vertex.cc[0] = (c[0].vertex.cc[0] + c[1].vertex.cc[0])/2;

                c[v1].vertex.cc[0] = (c[3].vertex.cc[0] + c[2].vertex.cc[0])/2;

                changed = 1;

            }

            while (c[v0].vertex.p.y > pfs->rect.q.y+MIDPOINT_ACCURACY && c[v1].vertex.p.y > pfs->rect.q.y+MIDPOINT_ACCURACY);

            if (0)

            {

check_rot_bottom_quarter:

                /* At this point we know that the condition for the following loop is true, so it

                 * can be a do...while rather than a while. */

                do

                {

                    /* Let's form (Y coords only):

                     *

                     * c[v0].vertex     c[v1].vertex

                     * x                x

                     * x                x

                     * m0               m1

                     * c[v0^1].vertex   c[v1^1].vertex

                     */

                    m0 = quarterpoint(c[v0].vertex.p.y, c[v0^1].vertex.p.y);

                    if (m0 <= pfs->rect.q.y+QUARTERPOINT_ACCURACY)

                        break;

                    m1 = quarterpoint(c[v1].vertex.p.y, c[v1^1].vertex.p.y);

                    if (m1 <= pfs->rect.q.y+QUARTERPOINT_ACCURACY)

                        break;

                    /* So, we can completely discard the left hand half of the patch. */

                    c[v0].vertex.p.x = quarterpoint(c[v0].vertex.p.x, c[v0^1].vertex.p.x);

                    c[v0].vertex.p.y = m0;

                    c[v1].vertex.p.x = quarterpoint(c[v1].vertex.p.x, c[v1^1].vertex.p.x);

                    c[v1].vertex.p.y = m1;

                    c[v0].vertex.cc[0] = (c[v0].vertex.cc[0] + 3*c[v0^1].vertex.cc[0])/4;

                    c[v1].vertex.cc[0] = (c[v1].vertex.cc[0] + 3*c[v1^1].vertex.cc[0])/4;

                    changed = 1;

                }

                while (c[v0].vertex.p.y > pfs->rect.q.y+QUARTERPOINT_ACCURACY && c[v1].vertex.p.y > pfs->rect.q.y+QUARTERPOINT_ACCURACY);

            }

        }

    } while (changed);

    c[0].vertex.cc[1] = c[1].vertex.cc[1] =

                        c[2].vertex.cc[1] =

                        c[3].vertex.cc[1] = 0;

    make_other_poles(c);

    return patch_fill(pfs, c, NULL, NULL);

}

#undef midpoint

#undef quarterpoint

#undef MIDPOINT_ACCURACY

#undef QUARTERPOINT_ACCURACY

#define f_fits_in_fixed(f) f_fits_in_bits(f, fixed_int_bits)

static int

A_fill_region_floats(patch_fill_state_t *pfs1, corners_and_curves *cc, int depth)

{

    corners_and_curves sub[4];

    int code;

    if (depth == 32)

        return gs_error_limitcheck;

    if (depth > 0 &&

        f_fits_in_fixed(cc->corners[0].x) &&

        f_fits_in_fixed(cc->corners[0].y) &&

        f_fits_in_fixed(cc->corners[1].x) &&

        f_fits_in_fixed(cc->corners[1].y) &&

        f_fits_in_fixed(cc->corners[2].x) &&

        f_fits_in_fixed(cc->corners[2].y) &&

        f_fits_in_fixed(cc->corners[3].x) &&

        f_fits_in_fixed(cc->corners[3].y))

    {

        cc->curve[0].vertex.p.x = float2fixed(cc->corners[0].x);

        cc->curve[0].vertex.p.y = float2fixed(cc->corners[0].y);

        cc->curve[1].vertex.p.x = float2fixed(cc->corners[1].x);

        cc->curve[1].vertex.p.y = float2fixed(cc->corners[1].y);

        cc->curve[2].vertex.p.x = float2fixed(cc->corners[2].x);

        cc->curve[2].vertex.p.y = float2fixed(cc->corners[2].y);

        cc->curve[3].vertex.p.x = float2fixed(cc->corners[3].x);

        cc->curve[3].vertex.p.y = float2fixed(cc->corners[3].y);

        return subdivide_patch_fill(pfs1, cc->curve);

    }

    /* We have patches with corners:

     *  0  1

     *  3  2

     * We subdivide these into 4 smaller patches:

     *

     *  0   10   1     Where 0123 are corners

     *   [0]  [1]      [0][1][2][3] are patches.

     *  3   23   2

     *  0   10   1

     *   [3]  [2]

     *  3   23   2

     */

    sub[0].corners[0].x = cc->corners[0].x;

    sub[0].corners[0].y = cc->corners[0].y;

    sub[1].corners[1].x = cc->corners[1].x;

    sub[1].corners[1].y = cc->corners[1].y;

    sub[2].corners[2].x = cc->corners[2].x;

    sub[2].corners[2].y = cc->corners[2].y;

    sub[3].corners[3].x = cc->corners[3].x;

    sub[3].corners[3].y = cc->corners[3].y;

    sub[1].corners[0].x = sub[0].corners[1].x = (cc->corners[0].x + cc->corners[1].x)/2;

    sub[1].corners[0].y = sub[0].corners[1].y = (cc->corners[0].y + cc->corners[1].y)/2;

    sub[3].corners[2].x = sub[2].corners[3].x = (cc->corners[2].x + cc->corners[3].x)/2;

    sub[3].corners[2].y = sub[2].corners[3].y = (cc->corners[2].y + cc->corners[3].y)/2;

    sub[3].corners[0].x = sub[0].corners[3].x = (cc->corners[0].x + cc->corners[3].x)/2;

    sub[3].corners[0].y = sub[0].corners[3].y = (cc->corners[0].y + cc->corners[3].y)/2;

    sub[2].corners[1].x = sub[1].corners[2].x = (cc->corners[1].x + cc->corners[2].x)/2;

    sub[2].corners[1].y = sub[1].corners[2].y = (cc->corners[1].y + cc->corners[2].y)/2;

    sub[0].corners[2].x = sub[1].corners[3].x =

                          sub[2].corners[0].x =

                          sub[3].corners[1].x = (sub[0].corners[3].x + sub[1].corners[2].x)/2;

    sub[0].corners[2].y = sub[1].corners[3].y =

                          sub[2].corners[0].y =

                          sub[3].corners[1].y = (sub[0].corners[3].y + sub[1].corners[2].y)/2;

    sub[0].curve[0].vertex.cc[0] = sub[0].curve[3].vertex.cc[0] =

                                   sub[3].curve[0].vertex.cc[0] =

                                   sub[3].curve[3].vertex.cc[0] = cc->curve[0].vertex.cc[0];

    sub[1].curve[1].vertex.cc[0] = sub[1].curve[2].vertex.cc[0] =

                                   sub[2].curve[1].vertex.cc[0] =

                                   sub[2].curve[2].vertex.cc[0] = cc->curve[1].vertex.cc[0];

    sub[0].curve[1].vertex.cc[0] = sub[0].curve[2].vertex.cc[0] =

                                   sub[1].curve[0].vertex.cc[0] =

                                   sub[1].curve[3].vertex.cc[0] =

                                   sub[2].curve[0].vertex.cc[0] =

                                   sub[2].curve[3].vertex.cc[0] =

                                   sub[3].curve[1].vertex.cc[0] =

                                   sub[3].curve[2].vertex.cc[0] = (cc->curve[0].vertex.cc[0] + cc->curve[1].vertex.cc[0])/2;

    depth++;

    if (not_clipped_away(sub[0].corners, &pfs1->rect)) {

        code = A_fill_region_floats(pfs1, &sub[0], depth);

        if (code < 0)

            return code;

    }

    if (not_clipped_away(sub[1].corners, &pfs1->rect)) {

        code = A_fill_region_floats(pfs1, &sub[1], depth);

        if (code < 0)

            return code;

    }

    if (not_clipped_away(sub[2].corners, &pfs1->rect)) {

        code = A_fill_region_floats(pfs1, &sub[2], depth);

        if (code < 0)

            return code;

    }

    if (not_clipped_away(sub[3].corners, &pfs1->rect)) {

        code = A_fill_region_floats(pfs1, &sub[3], depth);

        if (code < 0)

            return code;

    }

    return 0;

}

#define midpoint(a,b)      ((a+b)/2)

#define quarterpoint(a,b)  ((a+3*b)/4)

static int

subdivide_patch_fill_floats(patch_fill_state_t *pfs, corners_and_curves *cc)

{

    double m0, m1;

    int v0, v1;

    int changed;

    if (pfs->rect.p.x >= pfs->rect.q.x || pfs->rect.p.y >= pfs->rect.q.y)

        return 0;

    /* On entry we have a patch:

     *   c[0].vertex  c[1].vertex

     *

     *   c[3].vertex  c[2].vertex

     *

     * Only the corners are set. The control points are not!

     *

     * BUT... in terms of spacial coords, it might be different...

     * They might be flipped on X, Y or both, giving:

     *  01 or 10 or 32 or 23

     *  32    23    01    10