/*
___
\3 / Input: integer x >= 0
\/ x = y Output: integer y, such that y^3 <= x < (y+1)^3
Cube Roots are computed three to five times faster with the CR function,
compared to the Nth Root function.
For small values it's obvious to use: " y = (uint)Math.Pow(x, 1d / 3)) ".
It takes ~125 ns, but for example with x = 4503569204744003 (a 52 bits number),
it returns 165139, wrong, it should be 165140. So errors have to be corrected.
Another option: the "Integer Cube Root" algorithm from "Hacker's Delight" (see Refs).
For a uint it takes ~25 ns, not too bad for a C# version, for a ulong ~330 ns.
Up to ulong.MaxValue there are ~2.5 million perfect cubes, they can be generated like below.
3th powers:
|-----|------|------|------|------| dy=growth of y, ddy=growth of dy, ...
| x | y=x^3| dy | ddy | dddy |
|-----|------|------|------|------| y=x^3 , y'=3*x^2 , y''=6*x , y'''=6
| 0 | 0 | | 0 | |
| | | 1 | | 6 | y=x^n ,,,,,,,,,,,,,, y''''...=n! A000142
| 1 | 1 | | 6 | |
| | | 7 | | 6 | y: The cubes A000578
| 2 | 8 | | 12 | | dy: Central hexagonal numbers A003215
| | | 19 | | 6 | ddy: Multiples of 6 A008588
| 3 | 27 | | 18 | | dddy: The six sequence A010722
| | | 37 | | 6 | ddddy: The zero sequence A000004
| 4 | 64 | | 24 | |
|-----|------|------|------|------| y=x^4 : A000583 , A010863 , A101103 , A005914 , A005917
private static void abc()
{
uint n = 0, a = 0, b = 1, c = 6;
Console.WriteLine("// {0,2}{1,4}{2,4}{3,3}", n, a, b, c); // 0 0 1 6
while (n < 7) // 1 1 7 12
{ // 2 8 19 18
n++; // 3 27 37 24
a += b; // 4 64 61 30
b += c; // 5 125 91 36
c += 6; // 6 216 127 42
Console.WriteLine("// {0,2}{1,4}{2,4}{3,3}", n, a, b, c); // 7 343 169 48
}
Console.ReadLine();
}
4th powers:
private static void abcd()
{
uint n = 0;
uint a = 0, b = 1, c = 14, d = 36;
Console.WriteLine("// {0,2}{1,5}{2,5}{3,4}{4,4}", n, a, b, c, d); // 0 0 1 14 36
while (n < 8) // 1 1 15 50 60
{ // 2 16 65 110 84
n++; // 3 81 175 194 108
a += b; // 4 256 369 302 132
b += c; // 5 625 671 434 156
c += d; // 6 1296 1105 590 180
d += 24; // 7 2401 1695 770 204
Console.WriteLine("// {0,2}{1,5}{2,5}{3,4}{4,4}", n, a, b, c, d); // 8 4096 2465 974 228
}
Console.ReadLine();
}
The initialisation value's: "a = 0, b = 1, c = 14, d = 36, 24" : A019538
It's getting a bit off topic, hence another topic: Power Addition Only.
Above leads to cro12 / CR12, "Hacker's Delight" hardware algorithms, and
a new "cro12", while writing ... this. It's ~1 ns faster than the one used.
private static uint cro12new(uint x) // |------|----|
{ // | x | ns |
uint y = 0, a = 0, b = 1, c = 0; // |------|----|
while (a < x) // | 0 | 4 |
{ // | 1 | 4 |
y++; // | 7 | 5 |
b += c; // | 63 | 7 |
a += b; // | 255 | 10 |
c += 6; // | 1000 | 13 |
} // | 1023 | 15 |
if (a != x) y--; // | 4095 | 24 |
return y; // |------|----|
}
|----------------------------------|
| Times in ns (Athlon X4, XP, 2GB) |
|----------|---------|-------------|
| | | output type |
| function | input |------|------|
| | | uint | Xint |
|----------|---------|------|------|
| cro12 | < 4096 | <=25 | <60 |
| cro32 | <= ~0u | ~25 | ~80 |
| cro64 | <= ~0uL | ~330 | ~360 |
| | | | |
| CR12 | < 4096 | <=35 | <130 |
| CR32 | <= ~0u | ~35 | ~135 |
| CR64 | <= ~0uL | ~340 | ~460 |
|----------------------------------|
Abbreviations:
ref reference
CR Cube Root & Remainder CR(9)={2,1}
CRO Cube Root Only CRO(9)= 2
RND Random number
~ bitwise not operator, or depending
on context, approximately
~0u uint.MaxValue
~0uL ulong.MaxValue
ns nanosecond (10^-9 s)
us microsecond (10^-6 s)
ms millisecond (10^-3 s)
Times:
|------------------------------------------------------------------------|
| Athlon X4, XP, 2GB |
|--------|--------|--------||----------|---------||------------|---------|
| | ns | ns || | us || | ms |
| X | CRO(X) | CR(X) || bits RND | CR(RND) || bits RND | CR(RND) |
|--------|--------|--------||----------|---------||------------|---------|
| 0 | 35 | 97 || 65 | 4 || 100.000 | 13 |
| 1 | 36 | 98 || 100 | 8 || 200.000 | 31 |
| 100 | 42 | 103 || 200 | 14 || 500.000 | 92 |
| 1000 | 47 | 112 || 500 | 27 || 1.000.000 | 237 |
| 4095 | 59 | 124 || 1000 | 38 || 2.000.000 | 620 |
| 4096 | 66 | 129 || 2000 | 59 || 5.000.000 | 2330 |
| ~0u | 83 | 142 || 5000 | 134 || 10.000.000 | 6250 |
| ~0u+1 | 329 | 395 || 10000 | 341 || 20.000.000 | 16500 |
| ~0uL | 382 | 474 || 20000 | 1060 || | |
| ~0uL+1 | 3870 | 3820 || 50000 | 4360 || | |
| | | || 100000 | 12400 || | |
|--------|--------|--------||----------|---------||------------|---------|
Refs: Hacker's Delight
Henry S. Warren, Jr.
ISBN 0-201-91465-4
8th Printing February 2008
http://www.hackersdelight.org
(11-2 Integer Cube Root)
Modern Computer Arithmetic
Richard Brent and Paul Zimmermann
version 0.1, October 2006
http://www.loria.fr/~zimmerma/mca/mca-0.1.pdf
(1.5.2 k-th Root (Cube Root))
*/
using System;
using System.Diagnostics;
using System.Threading.Tasks;
using Xint = System.Numerics.BigInteger;
class Cube_Root
{
private static Xint[] CR(Xint D)
{
if (D <= ~0uL) return D < 4096 ? CR12((uint)D) : D <= ~0u ? CR32((uint)D) : CR64((ulong)D);
int n = bL(D) / 6, m = 2 * n;
Xint BA = D & ((Xint.One << m) - 1); // |-------D------|
D >>= m; // | |
Xint C = D & ((Xint.One << n) - 1); // |-D-|-C-|--BA--|
D >>= n;
Xint[] R = CR(D); // R[0] = CubeRoot
Xint R0R0 = SQ(R[0]);
Xint[] Q = DQR(C + (R[1] << n), 3 * R0R0); // Q[0] = Quotient
Xint Q0Q0 = SQ(Q[0]);
R[1] = BA + (Q[1] << m) - MTP(Q0Q0, Q[0] + ((3 * R[0]) << n));
if (R[1] < 0)
{
R0R0 <<= m;
R0R0 += Q0Q0 + (MTP(Q[0], R[0]) << (n + 1));
R[0] <<= n;
R[0] += Q[0];
R[1] += 1 + 3 * (R0R0 - R[0]);
R[0] -= 1;
while (R[1] < 0)
{
R0R0 -= 1 + 2 * R[0];
R[1] += 1 + 3 * (R0R0 - R[0]);
R[0] -= 1;
}
}
else
{
R[0] <<= n;
R[0] += Q[0];
}
return R;
}
private static Xint[] CR12(uint x)
{
uint y = 0, z = 0, r = 0, s = 0;
while (r < x)
{
s = r;
z += y * 6;
r += z + 1;
y += 1;
}
return r == x ?
new Xint[] { y, 0 } : // no need to use s, replace "x - s" by "x - (r + z + 1)",
new Xint[] { y - 1, x - s }; // but it's much slower, why?
}
private static Xint[] CR32(uint x)
{
uint y = 0, z = 0, b = 0;
int s = 30;
while (s >= 0)
{
y *= 2;
z *= 4;
b = 3 * y + 3 * z + 1 << s;
s -= 3;
if (x >= b)
{
x -= b;
z += 2 * y + 1;
y += 1;
}
}
return new Xint[] { y, x };
}
private static Xint[] CR64(ulong x)
{
uint y = 0;
ulong z = 0, b = 0, bs = 0;
int s = 63;
while (s >= 0)
{
y *= 2;
z *= 4;
b = 3 * y + 3 * z + 1;
bs = b << s;
if (x >= bs && b == bs >> s)
{
x -= bs;
z += 2 * y + 1;
y += 1;
}
s -= 3;
}
return new Xint[] { y, x };
}
private static Xint CRO(Xint X)
{
if (X <= ~0ul) return X < 4096 ? cro12((uint)X) : X <= ~0u ? cro32((uint)X) : cro64((ulong)X);
Xint[] R = CR(X);
return R[0];
}
private static uint cro12(uint x)
{
uint y = 0, z = 0, r = 0;
while (r < x)
{
z += y * 6;
r += z + 1;
y += 1; // y > 1625 causes overflow of r
}
return r == x ? y : y - 1;
}
private static uint cro32(uint x)
{
uint y = 0, z = 0, b = 0;
int s = 30;
while (s >= 0)
{
y *= 2;
z *= 4;
b = 3 * y + 3 * z + 1 << s;
s -= 3;
if (x >= b)
{
x -= b;
z += 2 * y + 1;
y += 1;
}
}
return y;
}
private static uint cro64(ulong x)
{
uint y = 0;
ulong z = 0, b = 0, bs = 0;
int s = 63;
while (s >= 0)
{
y *= 2;
z *= 4;
b = 3 * y + 3 * z + 1;
bs = b << s;
if (x >= bs && b == bs >> s)
{
x -= bs;
z += 2 * y + 1;
y += 1;
}
s -= 3;
}
return y;
}
// Use the faster versions from: Nth Root ////////////////////////////////////////////////////
private static Xint SQ(Xint U) { return U * U; } //
private static Xint MTP(Xint U, Xint V) { return U * V; } //
private static Xint[] DQR(Xint U, Xint V) //
{ //
Xint[] QR = new Xint[2]; //
QR[0] = Xint.DivRem(U, V, out QR[1]); //
return QR; //
} //
// Use the faster versions: http://www.bigintegers.blogspot.com/2013/07/nth-root-power.html //
private static Stopwatch sw = new Stopwatch();
static void Main()
{
for (int i = 1; i < 100; i++)
{
Xint X = RND(8 * i);
//X = Xint.Pow(X, 3);
//X--;
//X++;
sw.Start(); Xint[] R = CR(X); sw.Stop();
if (Xint.Pow(R[0], 3) > X) Console.WriteLine("WRONG1");
if (Xint.Pow(R[0] + 1, 3) <= X) Console.WriteLine("WRONG2");
if (Xint.Pow(R[0], 3) + R[1] != X) Console.WriteLine("WRONG3");
}
Console.WriteLine(sw.ElapsedMilliseconds + " ms");
Console.ReadLine();
}
private static int bL(Xint U)
{
byte[] bytes = (U.Sign * U).ToByteArray();
int i = bytes.Length - 1;
return i * 8 | bitLengthMostSignificantByte(bytes[i]);
}
private static int bitLengthMostSignificantByte(byte b)
{
return b < 08 ? b < 02 ? b < 01 ? 0 : 1 :
b < 04 ? 2 : 3 :
b < 32 ? b < 16 ? 4 : 5 :
b < 64 ? 6 : 7;
}
private static int seed;
private static Xint RND(int n)
{
if (n < 2) return n;
if (seed == int.MaxValue) seed = 0; else seed++;
Random rand = new Random(seed);
byte[] bytes = new byte[(n + 15) >> 3];
rand.NextBytes(bytes);
int i = bytes.Length - 1;
bytes[i] = 0;
n = (i << 3) - n;
i--;
bytes[i] >>= n;
bytes[i] |= (byte)(128 >> n);
return new Xint(bytes);
}
private static Xint[] CRold(Xint D)
{
if (D < 64)
{
int d = (int)D;
if (d >= 27) return new Xint[] { 3, d - 27 };
if (d >= 08) return new Xint[] { 2, d - 08 };
if (d >= 01) return new Xint[] { 1, d - 01 };
return new Xint[] { 0, 0 };
}
int n = bL(D) / 6;
Xint Mask = (Xint.One << n) - 1;
Xint A = D & Mask; D >>= n;
Xint B = D & Mask; D >>= n;
Xint C = D & Mask; D >>= n;
Xint[] R = CRold(D); // R[0] = CubeRoot, R[1] = Remainder
Xint R0R0 = SQ(R[0]);
Xint[] Q = DQR(C + (R[1] << n), 3 * R0R0); // Q[0] = Quotient, Q[1] = Remainder
Xint Q0Q0 = SQ(Q[0]);
R[1] = A + ((B + (Q[1] << n)) << n) - MTP(Q0Q0, Q[0] + ((3 * R[0]) << n));
if (R[1] < 0)
{
R0R0 <<= n * 2;
R0R0 += Q0Q0 + (MTP(Q[0], R[0]) << (n + 1));
R[0] <<= n;
R[0] += Q[0];
R[1] += 1 + 3 * (R0R0 - R[0]);
R[0] -= 1;
while (R[1] < 0)
{
R0R0 -= 1 + 2 * R[0];
R[1] += 1 + 3 * (R0R0 - R[0]);
R[0] -= 1;
}
}
else
{
R[0] <<= n;
R[0] += Q[0];
}
return R;
}
}
