Re: square number

/*
Figure out if a number is an integral square
Still not perfect, due to floating point imprecision
but it should work most of the time.
*/
#include <math.h>
double isroot(long x)
{
double d = x;
double root = sqrt(x);
double floor_of_root = floor(root);
double ceil_of_root = ceil(root);
if ((floor_of_root * floor_of_root) == x)
return floor_of_root;
if ((ceil_of_root * ceil_of_root) == x)
return ceil_of_root;
return -1;
}
#ifdef UNIT_TEST
#include <stdio.h>
char string[32767];
int main(void)
{
puts("Enter a number:");
fflush(stdout);
if (fgets(string, sizeof string, stdin) != NULL) {
long question = atol(string);
double answer = isroot(question);
if (answer < 0)
printf("%ld is not a square number\n", question);
else
printf("%.0f * %.0f = %ld\n", answer, answer, question);
}
return 0;
}
#endif


Here is the solution to a related problem:
Given a number, find what number raised to its own power is equal to
that number.
So given 4, the answer is 2 and given 27 the answer is 3. But what if
you are given 5 or 32.123?

#include <math.h>
#include <float.h>
/*****************************************************************************/
/* Crude approximation for the perfect root of a number. We want to be
able */
/* to calculate what number, when raised to its own power, will be
equal to */
/* the input number. For instance, if given 27, the answer will be 3,
since */
/* 3 to the 3rd power is 27. This function will provide an estimate
for a */
/* root solving function to begin iteration (Two or three correct
digits). */
/*****************************************************************************/
long double PerfectRootGuess( long double x )
{
long double Guess=0;
long double Test;

/* Min is (1/e)^(1/e). For smaller x, there are two distinct
answers. */
long double Min = 0.69220062755534635386542199718278976149L;

/* Choose Max such that Max^Max <= LDBL_MAX */
long double Max = 1546.0e0L;
int n;
for ( n = 0; n < 16; n++ )
{
Guess = ( Min + Max ) * 0.5e0L;
Test = powl( Guess, Guess );
if ( Test > x )
Max = Guess;
else
Min = Guess;
}
return Guess;
}
/********************************************************************/
/* This is Schroder's Method A( Lamda=0, Omega=4 ) for x^x - c = 0. */
/* Reference: "On Infinitely Many Algorithms For Solving Equations" */
/* by Ernst Schroder, Translated by G. W. Stewart. */
/* ---------------------------------------------------------------- */
/* Copyright (C) 1996 by Dann Corbit */
/********************************************************************/
long double PerfectRootImprove( long double c, long double x )
{
long double t1 = x * x;
long double t2 = t1 * x;
long double t3 = logl( x );
long double t4 = t3 * t3;
long double t5 = t4 * t3;
long double t6 = t2 * t5;
long double t7 = t1 * t1;
long double t8 = powl( t7, x );
long double t10 = t8 * t1;
long double t11 = t10 * t3;
long double t12 = powl( t1, 2.0 * x );
long double t13 = t12 * t7;
long double t14 = t4 * t4;
long double t15 = t13 * t14;
long double t16 = t13 * t4;
long double t18 = powl( x, x + 4.0 );
long double t20 = c * c;
long double t21 = t20 * c;
long double t22 = t18 * t4 * t21;
long double t23 = t12 * t2;
long double t24 = t23 * t4;
long double t25 = t7 * t8;
long double t26 = powl( t1, x );
long double t27 = t7 * t26;
long double t28 = t3 * t20;
long double t29 = t27 * t28;
long double t30 = t4 * t20;
long double t31 = t27 * t30;
long double t32 = t7 * t4;
long double t33 = powl( t2, x );
long double t34 = t33 * c;
long double t35 = t32 * t34;
long double t36 = t32 * t8;
long double t37 = t7 * t14;
long double t38 = t37 * t8;
long double t40 = t27 * t14 * t20;
long double t41 = 28.0 * t6 * t8 + 8.0 * t11 - 14.0 * t15 - 84.0 *
t16
- 6.0 * t22 - 108.0 * t24 + 13.0 * t25 - 44.0 *
t29
- 66.0 * t31 - 66.0 * t35 + 78.0 * t36 + 13.0 * t38
- 11.0 * t40 - 14.0 * t13;
long double t42 = t7 * t5;
long double t43 = t42 * t8;
long double t44 = t12 * t1;
long double t45 = t7 * t3;
long double t46 = t45 * t34;
long double t47 = t2 * t26;
long double t48 = t47 * t30;
long double t50 = powl( x, x + 3.0 );
long double t52 = t50 * t4 * t21;
long double t54 = t50 * t3 * t21;
long double t56 = t2 * t3 * t34;
long double t58 = t27 * t5 * t20;
long double t60 = powl( x, x + 1.0 );
long double t61 = t60 * t21;
long double t64 = t26 * t3 * t1 * t20;
long double t65 = t42 * t34;
long double t66 = t27 * t20;
long double t68 = t26 * t1 * t20;
long double t69 = - 60.0 * t23 + 52.0 * t43 - 22.0 * t44 - 44.0 *
t46
- 54.0 * t48 + 11.0 * t10 - 18.0 * t52 - 24.0 *
t54
+ 72.0 * t56 - 44.0 * t58 + 2.0 * t61 + 8.0 * t64
- 44.0 * t65 - 11.0 * t66 + 11.0 * t68;
long double t71 = t8 * t2;
long double t72 = t71 * t3;
long double t74 = t2 * t4 * t34;
long double t75 = t47 * t28;
long double t76 = t37 * t34;
long double t78 = t7 * t33 * c;
long double t79 = t8 * x;
long double t80 = t33 * t1;
long double t82 = t80 * t3 * c;
long double t83 = t50 * t21;
long double t84 = t33 * t2;
long double t88 = powl( x, x + 2.0 );
long double t89 = t88 * t21;
long double t93 = t18 * t3 * t21;
long double t94 = t18 * t21;
long double t96 = t18 * t14 * t21;
long double t97 = 168.0 * t72 + 54.0 * t74 - 72.0 * t75 - 11.0 * t76
- 11.0 * t78 - 2.0 * t79 + 8.0 * t82 - 10.0 * t83
+ 12.0 * t84 * t5 * c - 11.0 * t89 + 70.0 * t71
- 12.0 * t6 * t26 * t20 - 4.0 * t93 - t94 - t96;
long double t99 = t88 * t3 * t21;
long double t100 = t45 * t8;
long double t105 = t18 * t5 * t21;
long double t106 = t13 * t3;
long double t107 = t13 * t5;
long double t108 = t23 * t3;
long double t109 = t47 * t20;
long double t110 = t84 * c;
long double t112 = t12 * t3 * t1;
long double t114 = t33 * x * c;
long double t116 = t26 * x * t20;
long double t117 = t71 * t4;
long double t118 = t80 * c;
long double t119 = - 8.0 * t99 + 52.0 * t100 - 24.0 * t23 * t5
- 4.0 * t50 * t5 * t21 - 4.0 * t105 - 56.0 *
t106
- 56.0 * t107 - 144.0 * t108 - 30.0 * t109
+ 30.0 * t110 - 16.0 * t112 + 6.0 * t114 - 6.0 *
t116
+ 126.0 * t117 + 11.0 * t118;
long double t123 = - 4.0 * t11 - 14.0 * t15 - 84.0 * t16 - 6.0 *
t22
- 36.0 * t24 + 13.0 * t25 - 44.0 * t29 - 66.0 *
t31
- 66.0 * t35 + 78.0 * t36 + 13.0 * t38 - 11.0 *
t40
- 14.0 * t13;
long double t124 = - 36.0 * t23 + 52.0 * t43 + 2.0 * t44 - 44.0 *
t46
- 18.0 * t48 - t10 - 6.0 * t52 - 12.0 * t54
+ 36.0 * t56 - 44.0 * t58 - 2.0 * t61 - 4.0 *
t64
- 44.0 * t65 - 11.0 * t66;
long double t126 = - t68 + 84.0 * t72 + 18.0 * t74 - 36.0 * t75
- 11.0 * t76 - 11.0 * t78 + 2.0 * t79 - 4.0 *
t82
- 6.0 * t83 + t89 + 42.0 * t71 - 4.0 * t93 -
t94;
long double t127 = - t96 + 4.0 * t99 + 52.0 * t100 - 4.0 * t105
- 56.0 * t106 - 56.0 * t107 - 72.0 * t108 - 18.0
* t109
+ 18.0 * t110 + 8.0 * t112 - 6.0 * t114 + 6.0 *
t116
+ 42.0 * t117 - t118;
return x * ( t41 + t69 + t97 + t119 ) / ( t123 + t124 + t126 + t127
);
}
/********************************************************************/
/* This is Schroder's Method A(Lamda=1/2, Omega=1) for x^x - c = 0. */
/* Reference: "On Infinitely Many Algorithms For Solving Equations" */
/* by Ernst Schroder, Translated by G. W. Stewart. */
/* ---------------------------------------------------------------- */
/* Copyright (C) 1996 by Dann Corbit */
/********************************************************************/
long double PerfectRootImproveSlow( long double c, long double x )
{
long double Num; /* The numerator */
long double Den; /* The denominator */
const long double Term = powl( x, x );
const long double Lamda = 0.5e0L;
Num = x * ( Term - c );
Den = x * ( Term * ( logl( x ) + 1.0e0L ) ) - Lamda * Num;
x -= ( Num / Den );
return x;
}
long double PerfectRoot( long double c )
{
int i;
long double x = PerfectRootImproveSlow( c, PerfectRootGuess( c ) );
for ( i = 0; i < 20; i++ )
if ( fabsl( c - powl( x, x ) ) / c >= LDBL_EPSILON*2.0e0L )
{
if ( c < 2.0e37 ) /* Then we can use the faster method... */
x = PerfectRootImprove( c, x );
else /* Much slower, but less prone to overflow. */
x = PerfectRootImproveSlow( c, x );
}
else
break;
return x;
}
#ifdef TEST
#include <stdio.h>
int main()
{

long double c;
long double Guess=0;
for ( c = 1.; c < 1000; c += 3. )
{
Guess = PerfectRoot( c );
printf( "N=%.0Lf, Guess=%.20LeL, Error=%.20LeL\n",
c, Guess, fabsl( c - powl( Guess, Guess ) )/c );
}
for ( c = 1.; c < LDBL_MAX; c *= 2. )
{
Guess = PerfectRoot( c );
printf( "N=%Le, Guess=%.20LeL, Error=%.20LeL\n",
c, Guess, fabsl( c - powl( Guess, Guess ) )/c );
}
/* Huge: */
c = 1.0e100L;
Guess = PerfectRoot( c );
printf( "N=%Le, Guess=%.20LeL, Error=%.20LeL\n",
c, Guess, fabsl( c - powl( Guess, Guess ) )/c );
/* Tiny: */
c = 0.75L;
Guess = PerfectRoot( c );
printf( "N=%Le, Guess=%.20LeL, Error=%.20LeL\n",
c, Guess, fabsl( c - powl( Guess, Guess ) )/c );

/* Huge Error: */
c = LDBL_MAX;
Guess = PerfectRoot( c );
printf( "N=%Le, Guess=%.20LeL, Error=%.20LeL\n",
c, Guess, fabsl( c - powl( Guess, Guess ) )/c );
/* Tiny: */
c = 0.2L;
Guess = PerfectRoot( c );
printf( "N=%Le, Guess=%.20LeL, Error=%.20LeL\n",
c, Guess, fabsl( c - powl( Guess, Guess ) )/c );

return 0;
}
#endif

.

Relevant Pages

  • Re: square number
    ... Figure out if a number is an integral square ... /* Crude approximation for the perfect root of a number. ... Guess = PerfectRoot(c); ...
    (sci.math)
  • Re: Perfect square efficient algorithm
    ... Test = powl(Guess, Guess); ... const long double Term = powl ... Guess = PerfectRoot(c); ...
    (sci.math)
Loading