Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/ptwXY_binaryOperators.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 /*
  2 # <<BEGIN-copyright>>
  3 # <<END-copyright>>
  4 */
  5 
  6 #include <cmath>
  7 #include <float.h>
  8 
  9 #include "ptwXY.h"
 10 
 11 #if defined __cplusplus
 12 namespace GIDI {
 13 using namespace GIDI;
 14 #endif
 15 
 16 static double ptwXY_mod2( double v, double m, int pythonMod );
 17 static nfu_status ptwXY_mul2_s_ptwXY( ptwXYPoints *n, ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double x1, double y1, double x2, double y2, int level );
 18 static nfu_status ptwXY_div_s_ptwXY( ptwXYPoints *n, ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double x1, double y1, double x2, double y2, 
 19     int level, int isNAN1, int isNAN2 );
 20 static ptwXYPoints *ptwXY_div_ptwXY_forFlats( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int safeDivide );
 21 static nfu_status ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXYPoints *ptwXY1, double x, double *y );
 22 /*
 23 ************************************************************
 24 */
 25 nfu_status ptwXY_slopeOffset( ptwXYPoints *ptwXY, double slope, double offset ) { 
 26 
 27     int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
 28     ptwXYPoint *p;
 29     ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
 30 
 31     if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
 32 
 33     for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = slope * p->y + offset;
 34     for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = slope * o->point.y + offset; 
 35     return( ptwXY->status );
 36 }   
 37 /*
 38 ************************************************************
 39 */
 40 nfu_status ptwXY_add_double( ptwXYPoints *ptwXY, double value ) { return( ptwXY_slopeOffset( ptwXY, 1., value ) ); }
 41 nfu_status ptwXY_sub_doubleFrom( ptwXYPoints *ptwXY, double value ) { return( ptwXY_slopeOffset( ptwXY,  1., -value ) ); }
 42 nfu_status ptwXY_sub_fromDouble( ptwXYPoints *ptwXY, double value ) { return( ptwXY_slopeOffset( ptwXY, -1.,  value ) ); }
 43 nfu_status ptwXY_mul_double( ptwXYPoints *ptwXY, double value ) { return( ptwXY_slopeOffset( ptwXY, value, 0. ) ); }
 44 nfu_status ptwXY_div_doubleFrom( ptwXYPoints *ptwXY, double value ) { 
 45 
 46     if( value == 0. ) {
 47         ptwXY->status = nfu_divByZero; }
 48     else {
 49         ptwXY_slopeOffset( ptwXY, 1. / value, 0. ); 
 50     }
 51     return( ptwXY->status );
 52 }
 53 nfu_status ptwXY_div_fromDouble( ptwXYPoints *ptwXY, double value ) {
 54 /*
 55 *   This does not do any infilling and it should?????????
 56 */
 57 
 58     int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
 59     ptwXYPoint *p;
 60     ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
 61 
 62     if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
 63     if( ptwXY->interpolation == ptwXY_interpolationOther ) return( nfu_otherInterpolation );
 64 
 65     for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) if( p->y == 0. ) ptwXY->status = nfu_divByZero;
 66     for( o = overflowHeader->next; o != overflowHeader; o = o->next ) if( o->point.y == 0. ) ptwXY->status = nfu_divByZero;
 67     if( ptwXY->status != nfu_divByZero ) {
 68         for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = value / p->y;
 69         for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = value / o->point.y; 
 70     }
 71     return( ptwXY->status );
 72 }
 73 /*
 74 ************************************************************
 75 */
 76 nfu_status ptwXY_mod( ptwXYPoints *ptwXY, double m, int pythonMod ) { 
 77 
 78     int64_t i, nonOverflowLength = ptwXY_getNonOverflowLength( ptwXY );
 79     ptwXYPoint *p;
 80     ptwXYOverflowPoint *o, *overflowHeader = &(ptwXY->overflowHeader);
 81 
 82     if( ptwXY->status != nfu_Okay ) return( ptwXY->status );
 83     if( m == 0 ) return( ptwXY->status = nfu_divByZero );
 84 
 85     for( i = 0, p = ptwXY->points; i < nonOverflowLength; i++, p++ ) p->y = ptwXY_mod2( p->y, m, pythonMod );
 86     for( o = overflowHeader->next; o != overflowHeader; o = o->next ) o->point.y = ptwXY_mod2( o->point.y, m, pythonMod );
 87     return( ptwXY->status );
 88 }
 89 /*
 90 ************************************************************
 91 */
 92 static double ptwXY_mod2( double v, double m, int pythonMod ) {
 93 
 94     double r = std::fmod( std::fabs( v ), std::fabs( m ) );
 95 
 96     if( pythonMod ) {
 97         if( ( v * m ) < 0. ) r = std::fabs( m ) - std::fabs( r );
 98         if( m < 0. ) r *= -1.; }
 99     else {
100         if( v < 0. ) r *= -1.;
101     }
102 
103     return( r );
104 }
105 /*
106 ************************************************************
107 */
108 ptwXYPoints *ptwXY_binary_ptwXY( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double v1, double v2, double v1v2, nfu_status *status ) {
109 
110     int64_t i;
111     int unionOptions = ptwXY_union_fill | ptwXY_union_mergeClosePoints;
112     double y;
113     ptwXYPoints *n;
114     ptwXYPoint *p;
115 
116     *status = nfu_otherInterpolation;
117     if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
118     if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
119     if( ( *status = ptwXY_areDomainsMutual( ptwXY1, ptwXY2 ) ) != nfu_Okay ) return( NULL );
120     if( ( ptwXY1->interpolation == ptwXY_interpolationFlat ) || ( ptwXY2->interpolation == ptwXY_interpolationFlat ) ) {
121         *status = nfu_invalidInterpolation;
122         if( ( ptwXY1->interpolation != ptwXY2->interpolation ) ) return( NULL );
123     }
124     if( ( n = ptwXY_union( ptwXY1, ptwXY2, status, unionOptions ) ) != NULL ) {
125         for( i = 0, p = n->points; i < n->length; i++, p++ ) {
126             if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, p->x, &y ) ) != nfu_Okay ) goto Err;
127             p->y = v1 * p->y + v2 * y + v1v2 * y * p->y;
128         }
129     }
130     return( n );
131 Err:
132     if( n ) ptwXY_free( n );
133     return( NULL );
134 }
135 /*
136 ************************************************************
137 */
138 ptwXYPoints *ptwXY_add_ptwXY( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status ) {
139 
140     ptwXYPoints *sum;
141 
142     if( ptwXY1->length == 0 ) {
143         sum = ptwXY_clone( ptwXY2, status ); }
144     else if( ptwXY2->length == 0 ) {
145         sum = ptwXY_clone( ptwXY1, status ); }
146     else {
147         sum = ptwXY_binary_ptwXY( ptwXY1, ptwXY2, 1., 1., 0., status );
148     }
149     return( sum );
150 }
151 /*
152 ************************************************************
153 */
154 ptwXYPoints *ptwXY_sub_ptwXY( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status ) {
155 
156     ptwXYPoints *diff;
157 
158     if( ptwXY1->length == 0 ) {
159         diff = ptwXY_clone( ptwXY2, status );
160         if( ( *status = ptwXY_neg( diff ) ) != nfu_Okay ) diff = ptwXY_free( diff ); }
161     else if( ptwXY2->length == 0 ) {
162         diff = ptwXY_clone( ptwXY1, status ); }
163     else {
164         diff = ptwXY_binary_ptwXY( ptwXY1, ptwXY2, 1., -1., 0., status );
165     }
166     return( diff );
167 }
168 /*
169 ************************************************************
170 */
171 ptwXYPoints *ptwXY_mul_ptwXY( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status ) {
172 
173     ptwXYPoints *mul;
174 
175     if( ptwXY1->length == 0 ) {
176         mul = ptwXY_clone( ptwXY1, status ); }
177     else if( ptwXY2->length == 0 ) {
178         mul = ptwXY_clone( ptwXY2, status ); }
179     else {
180         mul = ptwXY_binary_ptwXY( ptwXY1, ptwXY2, 0., 0., 1., status );
181     }
182     return( mul );
183 }
184 /*
185 ************************************************************
186 */
187 ptwXYPoints *ptwXY_mul2_ptwXY( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status ) {
188 
189     int64_t i, length;
190     ptwXYPoints *n = NULL;
191     int found;
192     double x1, y1, x2, y2, u1, u2, v1, v2, xz1 = 0, xz2 = 0, x;
193 
194     *status = nfu_otherInterpolation;
195     if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
196     if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
197     if( ( n = ptwXY_mul_ptwXY( ptwXY1, ptwXY2, status ) ) == NULL ) return( n );
198     if( ptwXY1->interpolation == ptwXY_interpolationFlat ) return( n );
199     if( ptwXY2->interpolation == ptwXY_interpolationFlat ) return( n );
200     length = n->length - 1;
201     if( length > 0 ) {
202         x2 = n->points[length].x;
203         for( i = length - 1; i >= 0; i-- ) {             /* Find and add y zeros not currently in n's. */
204             x1 = n->points[i].x;
205             if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x1, &u1 ) ) != nfu_Okay ) goto Err;
206             if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x2, &u2 ) ) != nfu_Okay ) goto Err;
207             if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x1, &v1 ) ) != nfu_Okay ) goto Err;
208             if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x2, &v2 ) ) != nfu_Okay ) goto Err;
209             found = 0;
210             if( u1 * u2 < 0 ) {
211                 xz1 = ( u1 * x2 - u2 * x1 ) / ( u1 - u2 );
212                 if( ( *status = ptwXY_setValueAtX( n, xz1, 0. ) ) != nfu_Okay ) goto Err;
213                 found = 1;
214             }
215             if( v1 * v2 < 0 ) {
216                 xz2 = ( v1 * x2 - v2 * x1 ) / ( v1 - v2 );
217                 if( ( *status = ptwXY_setValueAtX( n, xz2, 0. ) ) != nfu_Okay ) goto Err;
218                 found += 1;
219             }
220             if( found > 1 ) {
221                 x = 0.5 * ( xz1 + xz2 );
222                 if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x, &u1 ) ) != nfu_Okay ) goto Err;
223                 if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x, &v1 ) ) != nfu_Okay ) goto Err;
224                 if( ( *status = ptwXY_setValueAtX( n, x, u1 * v1 ) ) != nfu_Okay ) goto Err;
225             }
226             x2 = x1;
227         }
228 
229         if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
230         length = n->length;
231         x2 = n->points[n->length-1].x;
232         y2 = n->points[n->length-1].y;
233         for( i = n->length - 2; i >= 0; i-- ) {             /* Make interpolation fit accuracy. Work backwards so new points will not mess up loop. */
234             x1 = n->points[i].x;
235             y1 = n->points[i].y;
236             if( ( *status = ptwXY_mul2_s_ptwXY( n, ptwXY1, ptwXY2, x1, y1, x2, y2, 0 ) ) != nfu_Okay ) goto Err;
237             x2 = x1;
238             y2 = y1;
239         }
240         ptwXY_update_biSectionMax( n, (double) length );
241     }
242     return( n );
243 
244 Err:
245     if( n ) ptwXY_free( n );
246     return( NULL );
247 }
248 /*
249 ************************************************************
250 */
251 static nfu_status ptwXY_mul2_s_ptwXY( ptwXYPoints *n, ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double x1, double y1, double x2, double y2, int level ) {
252 
253     nfu_status status;
254     double u1, u2, v1, v2, x, y, yp, dx, a1, a2;
255 
256     if( ( x2 - x1 ) < ClosestAllowXFactor * DBL_EPSILON * ( std::fabs( x1 ) + std::fabs( x2 ) ) ) return( nfu_Okay );
257     if( level >= n->biSectionMax ) return( nfu_Okay );
258     level++;
259     if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x1, &u1 ) ) != nfu_Okay ) return( status );
260     if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x2, &u2 ) ) != nfu_Okay ) return( status );
261     if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x1, &v1 ) ) != nfu_Okay ) return( status );
262     if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, x2, &v2 ) ) != nfu_Okay ) return( status );
263     if( ( u1 == u2 ) || ( v1 == v2 ) ) return( nfu_Okay );
264     a1 = u1 * v1;
265     if( y1 == 0 ) a1 = 0.;                                  /* Fix rounding problem. */
266     a2 = u2 * v2;
267     if( y2 == 0 ) a2 = 0.;                                  /* Fix rounding problem. */
268     if( ( a1 == 0. ) || ( a2 == 0. ) ) {                    /* Handle special case of 0 where accuracy can never be met. */
269         x = 0.5 * ( x1 + x2 ); }
270     else {
271         if( ( a1 * a2 < 0. ) ) return( nfu_Okay );  /* Assume rounding error and no point needed as zero crossings are set in ptwXY_mul2_ptwXY. */
272         a1 = std::sqrt( std::fabs( a1 ) );
273         a2 = std::sqrt( std::fabs( a2 ) );
274         x = ( a2 * x1 + a1 * x2 ) / ( a2 + a1 );
275     }
276     dx = x2 - x1;
277     yp = ( u1 * v1 * ( x2 - x ) + u2 * v2 * ( x - x1 ) ) / dx;
278     y = ( u1 * ( x2 - x ) + u2 * ( x - x1 ) ) * ( v1 * ( x2 - x ) + v2 * ( x - x1 ) ) / ( dx * dx );
279     if( std::fabs( y - yp ) < std::fabs( y * n->accuracy ) ) return( nfu_Okay );
280     if( ( status = ptwXY_setValueAtX( n, x, y ) ) != nfu_Okay ) return( status );
281     if( ( status = ptwXY_mul2_s_ptwXY( n, ptwXY1, ptwXY2, x, y, x2, y2, level ) ) != nfu_Okay ) return( status );
282     status = ptwXY_mul2_s_ptwXY( n, ptwXY1, ptwXY2, x1, y1, x, y, level );
283     return( status );
284 }
285 /*
286 ************************************************************
287 */
288 ptwXYPoints *ptwXY_div_ptwXY( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int safeDivide ) {
289 
290     int isNAN1, isNAN2;
291     int64_t i, j, k, zeros = 0, length, iYs;
292     double x1, x2, y1, y2, u1, u2, v1, v2, y, xz, nan = nfu_getNAN( ), s1, s2;
293     ptwXYPoints *n = NULL;
294     ptwXYPoint *p;
295 
296     if( ( *status = ptwXY_simpleCoalescePoints( ptwXY1 ) ) != nfu_Okay ) return( NULL );
297     if( ( *status = ptwXY_simpleCoalescePoints( ptwXY2 ) ) != nfu_Okay ) return( NULL );
298     *status = nfu_otherInterpolation;
299     if( ptwXY1->interpolation == ptwXY_interpolationOther ) return( NULL );
300     if( ptwXY2->interpolation == ptwXY_interpolationOther ) return( NULL );
301     if( ( ptwXY1->interpolation == ptwXY_interpolationFlat ) || ( ptwXY2->interpolation == ptwXY_interpolationFlat ) )
302         return( ptwXY_div_ptwXY_forFlats( ptwXY1, ptwXY2, status, safeDivide ) );
303     
304     if( ( *status = ptwXY_areDomainsMutual( ptwXY1, ptwXY2 ) ) != nfu_Okay ) return( NULL );
305     if( ( n = ptwXY_union( ptwXY1, ptwXY2, status, ptwXY_union_fill | ptwXY_union_mergeClosePoints ) ) != NULL ) {
306         for( i = 0, p = n->points; i < n->length; i++, p++ ) {
307             if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, p->x, &y ) ) != nfu_Okay ) goto Err;
308             if( y == 0. ) {
309                 if( p->y == 0. ) {
310                     iYs = 0;
311                     y1 = 0.;
312                     y2 = 0.;
313                     if( i > 0 ) {
314                         if( ( *status = ptwXY_getSlopeAtX( ptwXY1, p->x, '-', &s1 ) ) != nfu_Okay ) {
315                             if( *status != nfu_XOutsideDomain ) goto Err;
316                             s1 = 0.;
317                         }
318                         if( ( *status = ptwXY_getSlopeAtX( ptwXY2, p->x, '-', &s2 ) ) != nfu_Okay ) goto Err;
319                         if( s2 == 0. ) {
320                             y1 = nan; }
321                         else {
322                             y1 = s1 / s2;
323                         } 
324                         iYs++;
325                     }
326                     if( i < ( n->length - 1 ) ) {
327                         if( ( *status = ptwXY_getSlopeAtX( ptwXY1, p->x, '+', &s1 ) ) != nfu_Okay ) {
328                             if( *status != nfu_XOutsideDomain ) goto Err;
329                             s1 = 0.;
330                         }
331                         if( ( *status = ptwXY_getSlopeAtX( ptwXY2, p->x, '+', &s2 ) ) != nfu_Okay ) goto Err;
332                         if( s2 == 0. ) {
333                             y2 = nan; }
334                         else {
335                             y2 = s1 / s2;
336                         }
337                         iYs++;
338                     }
339                     p->y = ( y1 + y2 ) / iYs; 
340                     if( nfu_isNAN( p->y ) ) zeros++; }
341                 else {
342                     if( !safeDivide ) {
343                         *status = nfu_divByZero;
344                         goto Err;
345                     }
346                     zeros++;
347                     p->y = nan;
348                 } }
349             else {
350                 p->y /= y;
351             }
352         }
353         length = n->length - 1;
354         if( length > 0 ) {
355             x2 = n->points[length].x;
356             for( i = length - 1; i >= 0; i-- ) {             /* Find and add y zeros and NAN not currently in n's. */
357                 x1 = n->points[i].x;
358                 if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x1, &u1 ) ) != nfu_Okay ) goto Err;
359                 if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x2, &u2 ) ) != nfu_Okay ) goto Err;
360                 if( ( *status = ptwXY_getValueAtX( ptwXY2, x1, &v1 ) ) != nfu_Okay ) goto Err;
361                 if( ( *status = ptwXY_getValueAtX( ptwXY2, x2, &v2 ) ) != nfu_Okay ) goto Err;
362                 if( u1 * u2 < 0 ) {
363                     xz = ( u1 * x2 - u2 * x1 ) / ( u1 - u2 );
364                     if( ( *status = ptwXY_setValueAtX( n, xz, 0. ) ) != nfu_Okay ) goto Err;
365                 }
366                 if( v1 * v2 < 0 ) {
367                     if( !safeDivide ) {
368                         *status = nfu_divByZero;
369                         goto Err;
370                     }
371                     zeros++;
372                     xz = ( v1 * x2 - v2 * x1 ) / ( v1 - v2 );
373                     if( ( *status = ptwXY_setValueAtX( n, xz, nan ) ) != nfu_Okay ) goto Err;
374                 }
375                 x2 = x1;
376             }
377             if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
378             length = n->length;
379             x2 = n->points[n->length-1].x;
380             y2 = n->points[n->length-1].y;
381             isNAN2 = nfu_isNAN( y2 );
382             for( i = n->length - 2; i >= 0; i-- ) {             /* Make interpolation fit accuracy. Work backwards so new points will not mess up loop. */
383                 x1 = n->points[i].x;
384                 y1 = n->points[i].y;
385                 isNAN1 = nfu_isNAN( y1 );
386                 if( !isNAN1 || !isNAN2 ) {
387                     if( ( *status = ptwXY_div_s_ptwXY( n, ptwXY1, ptwXY2, x1, y1, x2, y2, 0, isNAN1, isNAN2 ) ) != nfu_Okay ) goto Err;
388                 }
389                 x2 = x1;
390                 y2 = y1;
391                 isNAN2 = isNAN1;
392             }
393             ptwXY_update_biSectionMax( n, (double) length );
394             if( zeros ) {
395                 if( ( *status = ptwXY_simpleCoalescePoints( n ) ) != nfu_Okay ) goto Err;
396                 for( i = 0; i < n->length; i++ ) if( !nfu_isNAN( n->points[i].y ) ) break;
397                 if( nfu_isNAN( n->points[0].y ) ) {                     /* Special case for first point. */
398                     if( i == n->length ) {                              /* They are all nan's, what now? */
399                         zeros = 0;
400                         for( i = 0; i < n->length; i++ ) n->points[i].y = 0.; }
401                     else {
402                          n->points[0].y = 2. * n->points[i].y;
403                         zeros--;
404                     }
405                 }
406                 for( i = n->length - 1; i > 0; i-- ) if( !nfu_isNAN( n->points[i].y ) ) break;
407                 if( nfu_isNAN( n->points[n->length - 1].y ) ) {         /* Special case for last point. */
408                     n->points[n->length - 1].y = 2. * n->points[i].y;
409                     zeros--;
410                 }
411                 if( zeros ) {
412                     for( i = 0; i < n->length; i++ ) if( nfu_isNAN( n->points[i].y ) ) break;
413                     for( k = i + 1, j = i; k < n->length; k++ ) {
414                     if( nfu_isNAN( n->points[k].y ) ) continue;
415                         n->points[j] = n->points[k];
416                         j++;
417                     }
418                     n->length = j;
419                 }
420             }
421         }
422     }
423     return( n );
424 
425 Err:
426     if( n ) ptwXY_free( n );
427     return( NULL );
428 }
429 /*
430 ************************************************************
431 */
432 static nfu_status ptwXY_div_s_ptwXY( ptwXYPoints *n, ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, double x1, double y1, double x2, double y2, 
433         int level, int isNAN1, int isNAN2 ) {
434 
435     nfu_status status;
436     double u1, u2, v1, v2, v, x, y, yp, dx, a1, a2;
437 
438     if( ( x2 - x1 ) < ClosestAllowXFactor * DBL_EPSILON * ( std::fabs( x1 ) + std::fabs( x2 ) ) ) return( nfu_Okay );
439     if( level >= n->biSectionMax ) return( nfu_Okay );
440     level++;
441     if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x1, &u1 ) ) != nfu_Okay ) return( status );
442     if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x2, &u2 ) ) != nfu_Okay ) return( status );
443     if( ( status = ptwXY_getValueAtX( ptwXY2, x1, &v1 ) ) != nfu_Okay ) return( status );
444     if( ( status = ptwXY_getValueAtX( ptwXY2, x2, &v2 ) ) != nfu_Okay ) return( status );
445     if( isNAN1 ) {
446         x = 0.5 * ( x1 + x2 );
447         if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x, &u1 ) ) != nfu_Okay ) return( status );
448         if( ( status = ptwXY_getValueAtX( ptwXY2, x, &v1 ) ) != nfu_Okay ) return( status );
449         y = u1 / v1; }
450     else if( isNAN2 ) {
451         x = 0.5 * ( x1 + x2 );
452         if( ( status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY1, x, &u2 ) ) != nfu_Okay ) return( status );
453         if( ( status = ptwXY_getValueAtX( ptwXY2, x, &v2 ) ) != nfu_Okay ) return( status );
454         y = u2 / v2; }
455     else {
456         if( ( u1 == u2 ) || ( v1 == v2 ) ) return( nfu_Okay );
457         if( ( y1 == 0. ) || ( y2 == 0. ) ) {                    /* Handle special case of 0 where accuracy can never be met. */
458             x = 0.5 * ( x1 + x2 ); }
459         else {
460             if( ( u1 * u2 < 0. ) ) return( nfu_Okay );  /* Assume rounding error and no point needed. */
461             a1 = std::sqrt( std::fabs( u1 ) );
462             a2 = std::sqrt( std::fabs( u2 ) );
463             x = ( a2 * x1 + a1 * x2 ) / ( a2 + a1 );
464         }
465         dx = x2 - x1;
466         v = v1 * ( x2 - x ) + v2 * ( x - x1 );
467         if( ( v1 == 0. ) || ( v2 == 0. ) || ( v == 0. ) ) return( nfu_Okay );     /* Probably not correct, but I had to do something. */
468         yp = ( u1 / v1 * ( x2 - x ) + u2 / v2 * ( x - x1 ) ) / dx;
469         y = ( u1 * ( x2 - x ) + u2 * ( x - x1 ) ) / v;
470         if( std::fabs( y - yp ) < std::fabs( y * n->accuracy ) ) return( nfu_Okay );
471     }
472     if( ( status = ptwXY_setValueAtX( n, x, y ) ) != nfu_Okay ) return( status );
473     if( ( status = ptwXY_div_s_ptwXY( n, ptwXY1, ptwXY2, x, y, x2, y2, level, 0, isNAN2 ) ) != nfu_Okay ) return( status );
474     status = ptwXY_div_s_ptwXY( n, ptwXY1, ptwXY2, x1, y1, x, y, level, isNAN1, 0 );
475     return( status );
476 }
477 /*
478 ************************************************************
479 */
480 static ptwXYPoints *ptwXY_div_ptwXY_forFlats( ptwXYPoints *ptwXY1, ptwXYPoints *ptwXY2, nfu_status *status, int safeDivide ) {
481 
482     int64_t i;
483     ptwXYPoints *n = NULL;
484     ptwXYPoint *p;
485     double y;
486 
487     *status = nfu_invalidInterpolation;
488     if( ptwXY1->interpolation != ptwXY_interpolationFlat ) return( NULL );
489     if( ptwXY2->interpolation != ptwXY_interpolationFlat ) return( NULL );
490     if( ( n = ptwXY_union( ptwXY1, ptwXY2, status, ptwXY_union_fill | ptwXY_union_mergeClosePoints ) ) != NULL ) {
491         for( i = 0, p = n->points; i < n->length; i++, p++ ) {
492             if( ( *status = ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXY2, p->x, &y ) ) != nfu_Okay ) goto Err;
493             if( y == 0. ) {
494                 if( ( safeDivide ) && ( p->y == 0 ) ) {
495                     *status = nfu_divByZero;
496                     goto Err;
497                 } }
498             else {
499                 p->y /= y;
500             }
501         }
502     }
503     return( n );
504 
505 Err:
506     if( n ) ptwXY_free( n );
507     return( NULL );
508 }
509 /*
510 ************************************************************
511 */
512 static nfu_status ptwXY_getValueAtX_ignore_XOutsideDomainError( ptwXYPoints *ptwXY1, double x, double *y ) {
513 
514     nfu_status status = ptwXY_getValueAtX( ptwXY1, x, y );
515 
516     if( status == nfu_XOutsideDomain ) status = nfu_Okay;
517     return( status );
518 }
519 
520 #if defined __cplusplus
521 }
522 #endif
523