[PROJ] geodtest/planimeter test failure

[Charles Karney]

I recommend you start with the sin(-90) != -1 failure.  In particular
can you check the inputs

   x, (double)qd

and outputs

   r, q

of the call to remquo in sincosdx in geodesic.c line 179.  Buggy
implementations of remquo have been a problem in the past.

Thanks...

   --Charles

On 9/14/24 2:26 PM, Greg Troxel via PROJ wrote:
> I have several failing tests.  I'm going to post about them in small
> groups, because I suspect they might have a common cause.  I realize the
> problem could be in my toolchain.  I'll rerun paranoia and NetBSD's
> regression tests.
> 
> My first question is: Do regression tests all pass for everyone else, on
> all platforms they try to run them on?  If anyone has a clean pass on
> anything that isn't windows/mac/linux, please email me privately and
> I'll summarize.  If you have any failures, please post here.
> 
> My environment
>    NetBSD 10 (up to date with the branch)
>    gcc version 10.5.0 (nb2 20230710)
>    building with cmake in a build directory
>    9th gen i7, dell workstation (not an odd computer)
> 
> geodtest prints:
> 
>    checkEquals fails: 1.001875e+07 != 3.002268e+07 +/- 1
>    checkEquals fails: 1.001875e+07 != 2.002072e+07 +/- 1
>    Planimeter0 fail: 2
>    checkEquals fails: 7.612398e-07 != 1000000 +/- 0.01
>    Planimeter29 fail: 1
> 
> Looking at Planimeter0 and quoting sparsely, the failing tests are
> these, on the same inputs.
> 
>    double pd[3][2] = {{90, 0}, {0, 0}, {0, 90}};
> 
>    planimeter(&g, pd, 3, &perimeter, &area);
>    result += checkEquals(perimeter, 30022685, 1);
>    result += checkEquals(area, 63758202715511.0, 1);
> 
>    polylength(&g, pd, 3, &perimeter);
>    result += checkEquals(perimeter, 20020719, 1);
> 
> 
> Do others build with -ffloat-store to avoid extended precision?  (I
> didn't find any such instructions in the readme.)
> 
> Planimeter29:
> 
>    geod_init(&g, wgs84_a, wgs84_f);
>    geod_polygon_init(&p, 0);
>    geod_polygon_addpoint(&g, &p, 0, 0);
>    geod_polygon_addedge(&g, &p,  90, 1000);
>    geod_polygon_addedge(&g, &p,   0, 1000);
>    geod_polygon_addedge(&g, &p, -90, 1000);
>    geod_polygon_compute(&g, &p, 0, 1, &area, nullptr);
>    /* The area should be 1e6.  Prior to the fix it was 1e6 - A/2, where
>     * A = ellipsoid area. */
>    result += checkEquals(area, 1000000.0, 0.01);
> 
> 
> geodsigntest prints:
> 
>    Line 149: sin(-90) != -1 (1)
>    1 failure
> 
> due to
> 
>    #define checksincosd(x, s, c) do {            \
>        T sx, cx;                                 \
>        geod_sincosd(x, &sx, &cx);                \
>        if (equiv(s, sx)) {                       \
>          printf("Line %d: sin(%g) != %g (%g)\n", \
>                 __LINE__, x, s, sx);             \
>          ++n;                                    \
>        }                                         \
>        if (equiv(c, cx)) {                       \
>          printf("Line %d: cos(%g) != %g (%g)\n", \
>                 __LINE__, x, c, cx);             \
>          ++n;                                    \
>        }                                         \
>      } while (0)
> 
>       checksincosd(- 90.0, -1.0, +0.0);
> 
> 
> I'll keep looking, but any clues appreciated.
> _______________________________________________
> PROJ mailing list
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj

-- 
Charles Karney <karney at alum.mit.edu>
702 Prospect Ave
Princeton, NJ 08540-4037