/********************************************************************/
/* Schroder's Methods A and B for the calculation of cube roots.    */
/* Reference: "On Infinitely Many Algorithms For Solving Equations" */
/*            by Ernst Schroder, Translated by G. W. Stewart.       */
/* ---------------------------------------------------------------- */
/* Copyright (C) 1996 by Dann Corbit                                */
/********************************************************************/
#include "cbrts.h"
/********************************************************************/
/* This is Schroder's Method(A, Lamda=-1, Omega = 1) for x^3-k = 0  */
/********************************************************************/
long double CbrtAM1( long double x, long double c )
{
   long double x2 = x * x;
   return 3.0 * x2 * x2 / ( 4.0 * x2 * x + c );
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=0, Omega = 1) for x^3 - k = 0 */
/********************************************************************/
long double CbrtA01( long double x, long double c )
{
   long double x2 = x * x;
   return 3.0 * x2 * x2 / ( 4.0 * x2 * x + c );
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=1, Omega = 1) for x^3 - k = 0 */
/********************************************************************/
long double CbrtAP1( long double x, long double c )
{
   long double x2 = x * x;
   return ( 2.0 * x2 * x - c ) / x2 / 3;
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=-1, Omega = 1) for x^3-k = 0  */
/********************************************************************/
long double CbrtBM1( long double x, long double c )
{
   long double x3 = x * x * x;
   return x * ( x3 - 2.0 * c ) / ( 2.0 * x3 - c );
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=0, Omega = 1) for x^3 - k = 0 */
/********************************************************************/
long double CbrtB01( long double x, long double c )
{
   return  -3.0 * x * c / ( x * x * x - 2.0 * c );
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=1, Omega = 1) for x^3 - k = 0 */
/********************************************************************/
long double CbrtBP1( long double x, long double c )
{
   return x * ( 4.0 * c + x * x * x ) / c / 3;
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=-1, Omega = 2) for x^3-k = 0  */
/********************************************************************/
long double CbrtAM2( long double x, long double c )
{
   long double x2 = x * x;
   long double x3 = x2 * x;
   long double x4 = x2 * x2;
   long double c2 = c * c;
   return 3.0 * x4 * ( 2.0 * x3 - c ) / ( 10.0 * x4 * x2 + 2.0 * x3 * c + c2 );
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=0, Omega = 2) for x^3 - k = 0 */
/********************************************************************/
long double CbrtA02( long double x, long double c )
{
   long double x3 = x * x * x;
   return x * ( x3 - 2.0 * c ) / ( 2.0 * x3 - c );
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=1, Omega = 2) for x^3 - k = 0 */
/********************************************************************/
long double CbrtAP2( long double x, long double c )
{
   long double x2 = x * x;
   long double x3 = x2 * x;
   long double x4 = x2 * x2;
   long double c2 = c * c;
   return 1 / x2 * ( x4 * x2 - 7.0 * x3 * c + c2 ) / ( x3 - 2.0 * c ) / 3;
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=-1, Omega = 2) for x^3-k = 0  */
/********************************************************************/
long double CbrtBM2( long double x, long double c )
{
   long double x2 = x * x;
   long double x3 = x2 * x;
   long double x4 = x2 * x2;
   long double c2 = c * c;
   return 1 / x2 * ( x4 * x2 - 7.0 * x3 * c + c2 ) / ( x3 - 2.0 * c ) / 3;
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=0, Omega = 2) for x^3 - k = 0 */
/********************************************************************/
long double CbrtB02( long double x, long double c )
{
   long double x2 = x * x;
   long double x3 = x2 * x;
   long double x4 = x2 * x2;
   long double c2 = c * c;
   return  -3.0 * x * c * ( 2.0 * x3 - c ) / ( x4 * x2 - 7.0 * x3 * c + c2 );
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=1, Omega = 2) for x^3 - k = 0 */
/********************************************************************/
long double CbrtBP2( long double x, long double c )
{
   long double x2 = x * x;
   long double x3 = x2 * x;
   return x * ( x3 - 2.0 * c ) / ( 2.0 * x3 - c );
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=-1, Omega = 3) for x^3-k = 0  */
/********************************************************************/
long double CbrtAM3( long double x, long double c )
{
   long double x2 = x * x;
   long double x3 = x2 * x;
   long double x4 = x2 * x2;
   long double x6 = x4 * x2;
   long double x8 = x4 * x4;
   long double c2 = c * c;
   return x4 * ( 10.0 * x6 - 16.0 * x3 * c + c2 ) /
   ( 20.0 * x8 * x - 4.0 * x6 * c + 4.0 * x3 * c2 + c2 * c );
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=0, Omega = 3) for x^3 - k = 0 */
/********************************************************************/
long double CbrtA03( long double x, long double c )
{
   long double x2 = x * x;
   long double x4 = x2 * x2;
   long double x6 = x4 * x2;
   long double x3c = x2 * x * c;
   long double c2 = c * c;
   return x * ( 4.0 * x6 - 19.0 * x3c + 4.0 * c2 ) /
   ( 10.0 * x6 - 16.0 * x3c + c2 );
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=1, Omega = 3) for x^3 - k = 0 */
/********************************************************************/
long double CbrtAP3( long double x, long double c )
{
   long double x2 = x * x;
   long double x4 = x2 * x2;
   long double x6 = x4 * x2;
   long double c2 = c * c;
   long double x3c = x2 * x * c;
   return x * ( x6 - 16.0 * x3c + 10.0 * c2 ) /
   ( 4.0 * x6 - 19.0 * x3c + 4.0 * c2 );
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=-1, Omega = 3) for x^3-k = 0  */
/********************************************************************/
long double CbrtBM3( long double x, long double c )
{
   long double x2 = x * x;
   long double x4 = x2 * x2;
   long double x6 = x4 * x2;
   long double c2 = c * c;
   long double x3c = x2 * x * c;
   return x * ( x6 - 16.0 * x3c + 10.0 * c2 ) /
   ( -19.0 * x3c + 4.0 * c2 + 4.0 * x6 );
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=0, Omega = 3) for x^3 - k = 0 */
/********************************************************************/
long double CbrtB03( long double x, long double c )
{
   long double x2 = x * x;
   long double x4 = x2 * x2;
   long double x6 = x4 * x2;
   long double x3c = x2 * x * c;
   long double c2 = c * c;
   return  -1 / x2 * c * ( 10.0 * x6 - 16.0 * x3c + c2 ) /
   ( x6 - 16.0 * x3c + 10.0 * c2 );
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=1, Omega = 3) for x^3 - k = 0 */
/********************************************************************/
long double CbrtBP3( long double x, long double c )
{
   long double x2 = x * x;
   long double x4 = x2 * x2;
   long double x6 = x4 * x2;
   long double c2 = c * c;
   long double x3c = x2 * x * c;
   return x * ( 4.0 * x6 - 19.0 * x3c + 4.0 * c2 ) / ( 10.0 * x6 - 16.0 * x3c + c2 );
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=-1, Omega = 4) for x^3-k = 0  */
/********************************************************************/
long double CbrtAM4( long double x, long double c )
{
   long double x2 = x * x;
   long double x3 = x2 * x;
   long double x4 = x2 * x2;
   long double x6 = x4 * x2;
   long double x8 = x4 * x4;
   long double c2 = c * c;
   long double c4 = c2 * c2;
   return 3.0 * x4 * x3 * ( 5.0 * x6 - 17.0 * x3 * c + 5.0 * c2 ) /
   ( 15.0 * x6 * c2 - 35.0 * x8 * x * c + 5.0 * x3 * c2 * c + 35.0 * x8 * x4 + c4 );
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=0, Omega = 4) for x^3 - k = 0 */
/********************************************************************/
long double CbrtA04( long double x, long double c )
{
   long double x2 = x * x;
   long double x3 = x2 * x;
   long double x4 = x2 * x2;
   long double x6 = x4 * x2;
   long double x8 = x4 * x4;
   long double c2 = c * c;
   return 1 / x2 * ( 5.0 * x8 * x - 45.0 * x6 * c + 30.0 * x3 * c2 - c2 * c ) /
   ( 5.0 * x6 - 17.0 * x3 * c + 5.0 * c2 ) / 3;
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=1, Omega = 4) for x^3 - k = 0 */
/********************************************************************/
long double CbrtAP4( long double x, long double c )
{
   long double x2 = x * x;
   long double x4 = x2 * x2;
   long double x8 = x4 * x4;
   long double x9 = x8 * x;
   long double c2 = c * c;
   long double c3 = c2 * c;
   long double x3c2 = x2 * x * c2;
   long double x6c = x4 * x2 * c;
   return x * ( x9 - 30.0 * x6c + 45.0 * x3c2 - 5.0 * c3 ) /
   ( 30.0 * x3c2 - 45.0 * x6c - c3 + 5.0 * x9 );
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=-1, Omega = 4) for x^3-k = 0  */
/********************************************************************/
long double CbrtBM4( long double x, long double c )
{
   long double x2 = x * x;
   long double x4 = x2 * x2;
   long double x8 = x4 * x4;
   long double x9 = x8 * x;
   long double c2 = c * c;
   long double c3 = c2 * c;
   long double x3c2 = x2 * x * c2;
   long double x6c = x4 * x2 * c;
   return x * ( x9 - 30.0 * x6c + 45.0 * x3c2 - 5.0 * c3 ) /
   ( -45.0 * x6c + 30.0 * x3c2 - c3 + 5.0 * x9 );
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=0, Omega = 4) for x^3 - k = 0 */
/********************************************************************/
long double CbrtB04( long double x, long double c )
{
   long double x2 = x * x;
   long double x3 = x2 * x;
   long double x4 = x2 * x2;
   long double x6 = x4 * x2;
   long double x8 = x4 * x4;
   long double c2 = c * c;
   return  -3.0 * x * c * ( 5.0 * x6 - 17.0 * x3 * c + 5.0 * c2 ) /
   ( x8 * x - 30.0 * x6 * c + 45.0 * x3 * c2 - 5.0 * c2 * c );
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=1, Omega = 4) for x^3 - k = 0 */
/********************************************************************/
long double CbrtBP4( long double x, long double c )
{
   long double x2 = x * x;
   long double x3 = x2 * x;
   long double x4 = x2 * x2;
   long double x6 = x4 * x2;
   long double x8 = x4 * x4;
   long double c2 = c * c;
   return 1 / x2 * ( 30.0 * x3 * c2 - 45.0 * x6 * c - c2 * c + 5.0 * x8 * x ) /
   ( 5.0 * x6 - 17.0 * x3 * c + 5.0 * c2 ) / 3;
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=-1, Omega = 5) for x^3-k = 0  */
/********************************************************************/
long double CbrtAM5( long double x, long double c )
{
   long double x2 = x * x;
   long double x3 = x2 * x;
   long double x4 = x2 * x2;
   long double x6 = x4 * x2;
   long double x7 = x4 * x3;
   long double x8 = x4 * x4;
   long double x9 = x8 * x;
   long double c2 = c * c;
   long double c3 = c2 * c;
   long double c4 = c2 * c2;
   return 3.0 * x7 * ( 7.0 * x9 - 42.0 * x6 * c + 30.0 * x3 * c2 - 2.0 * c3 ) /
   ( -126.0 * x8 * x4 * c + 14.0 * x6 * c3 + 70.0 * x9 * c2 + 6.0 * c4 * x3 + 56.0 * x8 * x7 + c4 * c );
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=0, Omega = 5) for x^3 - k = 0 */
/********************************************************************/
long double CbrtA05( long double x, long double c )
{
   long double x2 = x * x;
   long double x4 = x2 * x2;
   long double x8 = x4 * x4;
   long double x9 = x8 * x;
   long double c2 = c * c;
   long double c3 = c2 * c;
   long double x3c2 = x2 * x * c2;
   long double x6c = x4 * x2 * c;
   return x * ( 2.0 * x9 - 30.0 * x6c + 42.0 * x3c2 - 7.0 * c3 ) / ( 7.0 * x9 - 42.0 * x6c + 30.0 * x3c2 - 2.0 * c3 );
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=1, Omega = 5) for x^3 - k = 0 */
/********************************************************************/
long double CbrtAP5( long double x, long double c )
{
   long double x2 = x * x;
   long double x3 = x2 * x;
   long double x4 = x2 * x2;
   long double x6 = x4 * x2;
   long double x8 = x4 * x4;
   long double x9 = x8 * x;
   long double c2 = c * c;
   long double c3 = c2 * c;
   long double c4 = c2 * c2;
   return 1 / x2 * ( x4 * x8 - 50.0 * x9 * c + 141.0 * x6 * c2 - 50.0 * x3 * c3 + c4 ) /
   ( -30.0 * x6 * c - 7.0 * c3 + 42.0 * x3 * c2 + 2.0 * x9 ) / 3;
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=-1, Omega = 5) for x^3-k = 0  */
/********************************************************************/
long double CbrtBM5( long double x, long double c )
{
   long double x2 = x * x;
   long double x3 = x2 * x;
   long double x4 = x2 * x2;
   long double x6 = x4 * x2;
   long double x8 = x4 * x4;
   long double x9 = x8 * x;
   long double c2 = c * c;
   long double c3 = c2 * c;
   long double c4 = c2 * c2;
   return 1 / x2 * ( x8 * x4 - 50.0 * x9 * c + 141.0 * x6 * c2 - 50.0 * x3 * c3 + c4 ) /
   ( -30.0 * x6 * c - 7.0 * c3 + 42.0 * x3 * c2 + 2.0 * x9 ) / 3;
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=0, Omega = 5) for x^3 - k = 0 */
/********************************************************************/
long double CbrtB05( long double x, long double c )
{
   long double x2 = x * x;
   long double x3 = x2 * x;
   long double x4 = x2 * x2;
   long double x6 = x4 * x2;
   long double x8 = x4 * x4;
   long double x9 = x8 * x;
   long double c2 = c * c;
   long double c3 = c2 * c;
   long double c4 = c2 * c2;
   return  -3.0 * x * c * ( 7.0 * x9 - 42.0 * x6 * c + 30.0 * x3 * c2 - 2.0 * c3 ) /
   ( x8 * x4 - 50.0 * x9 * c + 141.0 * x6 * c2 - 50.0 * x3 * c3 + c4 );
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=1, Omega = 5) for x^3 - k = 0 */
/********************************************************************/
long double CbrtBP5( long double x, long double c )
{
   long double x2 = x * x;
   long double x4 = x2 * x2;
   long double x8 = x4 * x4;
   long double x9 = x8 * x;
   long double c2 = c * c;
   long double c3 = c2 * c;
   long double x3c2 = x2 * x * c2;
   long double x6c = x4 * x2 * c;
   return x * ( 2.0 * x9 - 30.0 * x6c + 42.0 * x3c2 - 7.0 * c3 ) / ( 7.0 * x9 - 42.0 * x6c + 30.0 * x3c2 - 2.0 * c3 );
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=-1, Omega = 6) for x^3-k = 0  */
/********************************************************************/
long double CbrtAM6( long double x, long double c )
{
   long double x2 = x * x;
   long double x3 = x2 * x;
   long double x4 = x2 * x2;
   long double x6 = x4 * x2;
   long double x7 = x4 * x3;
   long double x8 = x4 * x4;
   long double x9 = x8 * x;
   long double x12 = x8 * x4;
   long double x16 = x8 * x8;
   long double c2 = c * c;
   long double c3 = c2 * c;
   long double c4 = c2 * c2;
   return x7 * ( 28.0 * x12 - 266.0 * x9 * c + 357.0 * x6 * c2 - 77.0 * x3 * c3 + c4 ) /
   ( 84.0 * x16 * x2 + c4 * c2 + 7.0 * c4 * c * x3 + 301.0 * x12 * c2 - 336.0 * x8 * x7 * c + 7.0 * x9 * c3 + 21.0 * c4 * x6 );
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=0, Omega = 6) for x^3 - k = 0 */
/********************************************************************/
long double CbrtA06( long double x, long double c )
{
   long double x2 = x * x;
   long double x4 = x2 * x2;
   long double x8 = x4 * x4;
   long double x12 = x8 * x4;
   long double c2 = c * c;
   long double c4 = c2 * c2;
   long double x9c = x8 * x * c;
   long double x6c2 = x4 * x2 * c2;
   long double x3c3 = x2 * x * c2 * c;
   return x * ( 7.0 * x12 - 161.0 * x9c + 393.0 * x6c2 - 161.0 * x3c3 + 7.0 * c4 ) /
   ( 28.0 * x12 - 266.0 * x9c + 357.0 * x6c2 - 77.0 * x3c3 + c4 );
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=1, Omega = 6) for x^3 - k = 0 */
/********************************************************************/
long double CbrtAP6( long double x, long double c )
{
   long double x2 = x * x;
   long double x4 = x2 * x2;
   long double x8 = x4 * x4;
   long double x12 = x8 * x4;
   long double c2 = c * c;
   long double c4 = c2 * c2;
   long double x6c2 = x4 * x2 * c2;
   long double x3c3 = x2 * x * c2 * c;
   long double x9c = x8 * x * c;
   return x * ( x12 - 77.0 * x9c + 357.0 * x6c2 - 266.0 * x3c3 + 28.0 * c4 ) /
   ( 7.0 * x12 + 393.0 * x6c2 - 161.0 * x9c - 161.0 * x3c3 + 7.0 * c4 );
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=-1, Omega = 6) for x^3-k = 0  */
/********************************************************************/
long double CbrtBM6( long double x, long double c )
{
   long double x2 = x * x;
   long double x4 = x2 * x2;
   long double x8 = x4 * x4;
   long double x12 = x8 * x4;
   long double x9c = x8 * x * c;
   long double c2 = c * c;
   long double c4 = c2 * c2;
   long double x6c2 = x4 * x2 * c2;
   long double x3c3 = x2 * x * c2 * c;
   return x * ( x12 - 77.0 * x9c + 357.0 * x6c2 - 266.0 * x3c3 + 28.0 * c4 ) /
   ( 7.0 * x12 + 393.0 * x6c2 - 161.0 * x9c - 161.0 * x3c3 + 7.0 * c4 );
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=0, Omega = 6) for x^3 - k = 0 */
/********************************************************************/
long double CbrtB06( long double x, long double c )
{
   long double x2 = x * x;
   long double x4 = x2 * x2;
   long double x8 = x4 * x4;
   long double x12 = x8 * x4;
   long double c2 = c * c;
   long double c4 = c2 * c2;
   long double x9c = x8 * x * c;
   long double x6c2 = x4 * x2 * c2;
   long double x3c3 = x2 * x * c2 * c;
   return  -1 / x2 * c * ( 28.0 * x12 - 266.0 * x9c + 357.0 * x6c2 - 77.0 * x3c3 + c4 ) /
   ( x12 - 77.0 * x9c + 357.0 * x6c2 - 266.0 * x3c3 + 28.0 * c4 );
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=1, Omega = 6) for x^3 - k = 0 */
/********************************************************************/
long double CbrtBP6( long double x, long double c )
{
   long double x2 = x * x;
   long double x4 = x2 * x2;
   long double x8 = x4 * x4;
   long double x12 = x8 * x4;
   long double c2 = c * c;
   long double c4 = c2 * c2;
   long double x9c = x8 * x * c;
   long double x3c3 = x2 * x * c2 * c;
   long double x6c2 = x4 * x2 * c2;
   return x * ( 7.0 * x12 + 393.0 * x6c2 - 161.0 * x9c - 161.0 * x3c3 + 7.0 * c4 ) /
   ( 28.0 * x12 - 266.0 * x9c + 357.0 * x6c2 - 77.0 * x3c3 + c4 );
}