//
// IntBrick64.c
//
// This program is very similar to the IntBrickDbl.c program except that 64
// bit long long variables are used for the p*q*k integer calculations. In the
// earlier IntBrickDbl program only 32 bit integers were used for these
// calculations, and 64 bit integers allow greater precision.
//
// As in the original IntBrick.c program, this program is dedicated to trying
// to find a solution to the integer brick problem. The problem is to find
// integer dimensions to a rectangular solid such that all diagonals (including
// the main diagonal thru the center) are integer numbers.
//
// The program uses the following algorithm: If a solution exists, then a
// primitive solution must have an odd integer dimension for one edge of the
// brick and even integers for the other two dimensions. In the program, this
// "odd" side is referenced by the variable "Xodd", and trials are generated
// for X = 1,3,5,7, etc. For each odd X, all integer right triangles that use
// "X" are generated, and the other side of each of these triangles is placed
// in a list of candidates for the other dimensions of the brick.
//
// Note: If a solution exists, then the list will also include the side of the
// right triangle thru the center of the brick. The side of this center
// triangle is also the diagonal of the brick face with even edges. Since this
// diagonal squared must equal the sum of the squares of its sides, and these
// other two edges are included in the generated list; we just search the list
// to see if the sum of the squares of any two elements is equal to the square
// of some third element.
//
// To generate all possible right triangles that use a given odd number as one
// side, the following algebraic relation is used:
// X,Y = sides Z = hypotenuse
// = the 3 sides of a right triangle
//
// X = (P-Q)(P+Q)(K) (Current trial - odd number)
// Y = 2(P)(Q)(K) (Trial other edge from pythagorean theorem)
// Z = (P^2 + Q^2)(K) (Not used by program)
// (P, Q, K = integers)
//
// For each trial "odd X", X is partitioned into all possible three divisors
// (D1*D2*D3 = X) such that D1 = P-Q , D2 = P+Q, and D3 = K. Any permutation
// of D1,D2,D3 such that the second term is larger than the first will
// generate a valid right triangle that is then added to the list of edge
// candidates for the brick.
//
// All possible D1,D2,D3 combinations that evenly divide "X" are considered.
// When the list is complete, a routine is called to search the list. If the
// sum of the squares of any two list entries is equal to the square of some
// other entry, then a solution would be found.
//
// Program code uses several techniques to speed execution. While each trial
// "X" could be partitioned by dividing by all odd numbers <= the sqrt. of
// "X", this process is slow when "X" is large; hence divisors are found using
// a sieve algorithm and stored in a very large sparse array.
//
// 1) Generate all divisors for: 0 < odd Xs < BlockSize
// then process each odd X in this block.
// 2) Generate all divisors for: BlockSize < odd Xs < 2 x BlockSize
// then process each odd X in this block.
// 3) etc.
//
// As divisors are found they are stored as nodes in an array that would be
// "RowLimit" x "XoddPerBlock" if declared explicitly. Each col. in this sparse
// array is a link list of divisors (candidates for D1, D2, D3). Each odd X
// has its own link list starting with ColHeaders[].
//
// For each odd X, as right triangles are found, each edge (and the 64 low
// bits of edge^2) is stored in arrays for subsequent lookup/matching. Hashing
// is used to speed the lookups.
//
// A search segment will end at 8.001E9 above the starting point. The 8.001E9
// limit results from the usage of unsigned integers in the FirstCol[] array.
// Used for generating the 4.0E9 columns of odd numbers for the divisor nodes.
//
// Search is complete out to one trillion (1.0E12). If you want to extend the
// search beyond 1.0E12, "#define RowLimit" and the Down[] and Divisor[] arrays
// will have to be extended. Use the status report (output by the program) as a
// guide for other changes to the code.
//
//*****************************************************************************
#include <stdheaders.h> // The usual stdio.h, etc.
#define RowLimit 500032 // This will allow searches for "odd X" up
// to (2 * RowLimit)^2 = 1 Trillion.
// (= 1.0E+12 ).
#define XoddPerBlock 50000 // Can change this to alter memory used.
#define BlockSize 100000.0 // If "XoddPerBlock" is changed, also
// change "BlockSize" to 2 times the above,
// but floating point. Note: If BlockSize
// is > 300,000, then "LastRow" calc in
// InitSys() will have to be changed.
int main(void); // Define functions. Descriptions of what
void InitSys(void); // they do are included in the functions.
void MakeDivs(void);
void MakeTriangles(void);
void AddEdge(double, double, double);
void SearchEdges(void);
// All arrays are global. These 4 arrays
// are used by the sieve algorithm to
// generate the divisors for each odd X.
unsigned FirstCol[RowLimit]; // Each member = the 1st col used in a row
unsigned ColHeads[XoddPerBlock]; // Link list heads for "Xodd" divisors
unsigned Down[332000]; // Ptrs to next node down for an Xodd list
double Divisor[332000]; // Divisor for a node (D1,D2,etc.)
// For each trial "odd edge", the divisors
// are copied from the above to here.
double DivsList[2000] = {1.0}; // List of divisors for an "Xodd".
// DivsList[0] must always be 1.0.
int MaxDivNodes; // Used for range checking.
// Each of the 2nd sides for all possible right triangles is a candidate
// for the other two dimensions of the brick and the diagonal on the
// even surface of the brick. Information about them is stored in the
// following arrays.
// If the sum of the squares of any two edges (array members) is equal to
// the square of some other member, then the brick problem is solved. For
// all possible sums, a hashing system is used to speed up the search for
// this 3rd member.
double YforTrials[100000]; // Candidates for edges.
long long unsigned YsqrLow[100000]; // 64 low bits of the above squared. Does
// not include the 4 trailing "0" bits.
unsigned HashLink[100000]; // L.Lists for the above that share the same hash
unsigned HashHead[16384]; // Start of above linked lists. Hash index is bits
// 18 to 31 of YsqrLow.
char Databuff[80]; // Buffer for user input.
// Other global variables
double Xbase; // Increments at "BlockSize" per iteration up to
// (2*RowLimit)^2
double Xodd; // Odd side for brick - Try 1,3,5,7, etc.
double MaxD1; // Upper limit for 1st divisor ( < cbrt(Xodd) )
unsigned NbrOfDiv; // Nbr. of divisors for current "Xodd"
int MaxRow; // Index of last member in FirstCol[] that
// will be used for current odd X.
int LastRow; // Used to find MaxRow
int Div16, NotDiv16; // Ptrs. into "Y" candidate list ( YforTrials[] )
// Div16 indexes "Y" edges divisible by 16
// and works down from 24999.
// NotDiv16 indexes "Y" edges that are not divs.
// by 16 and works up from 25000.
int MaxNbrDiv = 1; // Largest nbr. of divisors so far.
int MinDiv16 = 25000; // Lowest and highest indexes so far
int MaxNot16 = 24999; // in the "Y" candidate lists.
double XbaseStart; // User can choose where to start the search
double XbaseEnd; // Will stop at XbaseStart + 8,001,000,000
int ExtCount = 0; // Count external solutions
int main(void) { // Find bricks with integer sides and diagonals
int i, j, col; // Misc. loop controls
double ResetD1atX; // Always == (MaxD1+2)^2 * (MaxD1+4)
InitSys(); // Initialize FirstCol[] & MaxD1
ResetD1atX = (MaxD1 + 2.0)*(MaxD1 + 2.0)*(MaxD1 + 4.0);
// Start main loop
for (Xbase = XbaseStart; Xbase < XbaseEnd; Xbase += BlockSize) {
if (fmod(Xbase, 1000000.0) == 0.0) { // Status report
printf("\nWorking on Xbase = %'.0f\n", Xbase);
printf("LastRow is up to %'d (Must stay < %'d)\n", LastRow, RowLimit);
printf("MaxNbrDiv = %d (Must stay < 2,000)\n", MaxNbrDiv);
printf("MaxDivNodes = %'d (Must stay under 332,000)\n", MaxDivNodes);
printf("MinDiv16 = %'d (Must stay > 0)\n", MinDiv16);
printf("MaxNot16 = %'d (Must stay < 100,000)\n", MaxNot16);
}
MakeDivs(); // Generate all the divisors for each of the "XoddPerBlock"
// odd numbers in the current block (BlockSize).
// For each odd "X" in the block
for (col = 0; col < XoddPerBlock; col++) {
Xodd += 2.0; // Next trial odd edge = 1,3,5,etc.
// Xodd = Xbase + 2.0 * col + 1.0
if (Xodd >= ResetD1atX) { // Test to see if limit for
MaxD1 += 2.0; // 1st divisor needs update.
ResetD1atX = (MaxD1 + 2.0)*(MaxD1 + 2.0)*(MaxD1 + 4.0);
}
j = 1; // Count divisors for current odd X
// Get list of divisors of X
for (i = ColHeads[col]; i; i = Down[i])
DivsList[j++] = Divisor[i]; // Note: DivsList[0] always = 1.0
if (j < 2) // Need at least 2 for sol.
continue;
if (j > MaxNbrDiv) // Range checking
MaxNbrDiv = j;
DivsList[j] = Xodd; // End of Divisors.
NbrOfDiv = j; // Save global var.
Div16 = 25000; // Indexes Y's divisible by 16
NotDiv16 = 24999; // Indexes Y's not divisible by 16
// Used to speed up search process
// Use the list of divisors of "X" to
MakeTriangles(); // generate all possible right
// triangles.
if (Div16 < MinDiv16) // Range checking
printf("At Xodd = %'.0f MinDiv16 = %'d (Must stay > 0)\n",
Xodd, MinDiv16 = Div16);
if (NotDiv16 > MaxNot16)
printf("At Xodd = %'.0f MaxNot16 = %'d (Must stay < 100,000)\n",
Xodd, MaxNot16 = NotDiv16);
// Then search these combinations
SearchEdges(); // to see if a solution exists.
} // End of cols. (odd Xs)
for (i = RowLimit - 1; i >= 0; i--) // Update 1st col ptrs for
FirstCol[i] -= XoddPerBlock; // next block of X's.
} // Repeat for next Xbase
printf("\nProgram ended normally at XbaseEnd = %'.0f\n", XbaseEnd);
return(0);
} // End program
//*****************************************************************************
//
// InitSys
//
// This routine asks the user for the "Xodd" starting base and then
// initializes the FirstCol[] array. For each row, (The actual trial divisor
// is 2 * row + 1) FirstCol[] tells where to place the first node in the
// sparse array for divisors of a number.
//
//*****************************************************************************
void InitSys(void) {
int i;
double UpperLimit;
double Idbl, gap;
double InitCol;
double Cols2skip, LastCol;
UpperLimit = 4.0 * RowLimit * RowLimit;
puts("Computer program by Bill Butler\n");
puts("You may start the search for integer bricks from any Xodd base");
printf("up to %'.0f. The number you enter will be rounded down to the\n",
UpperLimit);
printf("first %'.0f below your entry.).\n", BlockSize);
puts("The program will end at 8.001 billion above your start point.");
puts("\nEnter a number for \"Xbase start\".");
gets(Databuff);
XbaseStart = atof(Databuff);
XbaseStart -= fmod(XbaseStart, BlockSize);
XbaseEnd = XbaseStart + 8.001e+9;
Xodd = XbaseStart - 1.0; // Initialize the odd side of the brick.
// Will count up by 2s from here.
if (XbaseStart < 1.5e7) // LastRow is used to find the max row
LastRow = 2000; // ("MaxRow") that will be processed in
else // the current Xbase iteration. If Xbase
LastRow = sqrt(XbaseStart/4.0) + 20.0;
// is < 15 million, "MaxRow" will increase
// faster than the 20 increment in
// MakeDivisors(). Thus give it a
// headstart for low Xbase numbers.
Cols2skip = XbaseStart/2.0; // Will Start to the right of this many
// cols
LastCol = Cols2skip - 1.0; // Will use this if to right of other hits
// Calculations for the sieve algorithm.
// The content of FirstCol[] will be the calculated position for the following:
// Given: Candidates for "odd side" of brick = XbaseStart + 2*col + 1
// Given: A list of possible divisors = 2 * (index for FirstCol[]) + 3
// If the divisor is >= the sqrt of "odd side", then the
// first hit will be at the col for divisor^2. This will be at
// InitCol - Cols2skip (but use a limit to prevent numeric overflow)
// Else find the first number greater than "odd side" that generates an odd
// integer result when divided by the "divisor". The "Floor"
// calculation finds where this number will hit.
for (i = 0; i < RowLimit; i++) { // Set up 1st col for all divisors
Idbl = i; // Get floating point for calc.
InitCol = 2.0*Idbl*Idbl + 6.0*Idbl + 4.0; // used for all divisors
if (InitCol > Cols2skip) // Limit is used to prevent unsigned
FirstCol[i] = fmin((InitCol - Cols2skip), 4294967000.0); // int overflow
else {
gap = 2.0 * Idbl + 3.0;
FirstCol[i] = gap * floor((LastCol + gap - InitCol)/gap) + InitCol - Cols2skip;
} // The above calcs check OK as shown in
} // "IntBrickFirstCol" Excel spreadsheet
// Calculate the upper limit
// for the first divisor
MaxD1 = floor(cbrt(XbaseStart)); // Get cube root for initial MaxD1
if (fmod(MaxD1, 2.0) == 0.0) // Make it an odd number.
MaxD1 -= 1.0;
}
//*****************************************************************************
//
// Make Divisors
//
// This routine generates all the divisors for the "XoddPerBlock" odd numbers
// in the current block starting at the current "Xbase". The FirstCol[] array
// contains a number giving the first col (corresponds to an odd X) that is
// evenly divisible by the divisor for the row (Divisor equals 2*row + 3). The
// next number that is evenly divisible is "gap" cols to the right where "gap"
// is equal to the divisor. At each location, a node is created to save the
// divisor. For each of the XoddPerBlock trial odd numbers in the current
// block, a link list of trial divisors will be generated.
//
//*****************************************************************************
void MakeDivs(void) {
int i, NodeIndx;
unsigned col, gap;
for (i = XoddPerBlock - 1; i >= 0; i--)
ColHeads[i] = 0; // Initialize col heads
// Skip rows that won't be used.
for (i = LastRow; FirstCol[i] > XoddPerBlock; i--);
MaxRow = i;
if ((LastRow - i) < 20) // Update starting point
LastRow = i + 20;
NodeIndx = 0; // Counts nodes in sparse matrix.
for (i = MaxRow; i >= 0; i--) { // Start with last row
gap = i; // Gap is equal to divisor for the
gap <<= 1; // row.
gap += 3;
// Each node is a divisor for the
// particular column (= an odd X)
for (col = FirstCol[i]; col < XoddPerBlock; col += gap) {
NodeIndx++; // Increment number of nodes
Divisor[NodeIndx] = gap; // Will convert to double;
Down[NodeIndx] = ColHeads[col]; // Add to link list
ColHeads[col] = NodeIndx;
} // End of nodes for this row
FirstCol[i] = col; // Save col for next Xbase
} // End of rows for this Xbase
if (NodeIndx > MaxDivNodes) // Note: The value of MaxDivNodes
MaxDivNodes = NodeIndx; // is output in the status report.
} // End of make divisors routine
//*****************************************************************************
//
// Make Triangles
//
// Given a list of divisors of X (List is in DivsList[]), this routine finds
// divisor triplets (D1 * D2 * D3 = odd X), and then passes valid permutations
// of each combination to routine AddEdge() where they are posted to the
// search arrays. Any permutation such that the second number "D2" is larger
// than the first "D1" will generate a valid right triangle.
//
//*****************************************************************************
void MakeTriangles(void) {
double D1, D2, D3, remainder;
int i, j;
// Loop control also gets 1st divs.
for (i = 0; (D1 = DivsList[i]) <= MaxD1; i++) {
remainder = Xodd/D1; // D2 * D3 must equal this
// Get divisors 2 and 3. Process
// loop for all divisors "D2" that
j = i; // are <= the sqrt(remainder)
D3 = remainder/(D2=DivsList[j]);
while (D2 <= D3) {
if (D3 == floor(D3)) { // If division was exact,
AddEdge(D1,D3,D2); // then add 1st combination
if ( (D2>D1) && (D3>D2) ) { // to list. Also
AddEdge(D1,D2,D3); // add others if
AddEdge(D2,D3,D1); // permutations are OK.
}
}
D3 = remainder/(D2 = DivsList[++j]);
} // End while
} // End "for i"
} // End MakeTriangles
//*****************************************************************************
//
// Add Edge
//
// The 3 numbers passed to this routine are known to be valid divisors for the
// odd side of the brick. ( D1 * D2 * D3 = Xodd) The routine breaks them down
// into components "p", "q", and "k"
// (p - q) * (p + q) * k = Xodd
// which are used to compute the "Y" side of a right triangle; which is thus a
// candidate for an even edge of an integer brick.
// 2 * p * q * k = "Y" edge.
// Each of of these edge candidates is added to a list for later processing.
//
// Data in each node in the lists consists of trial edge "Y" which is stored
// as a "double" variable, and the 64 low bits of Y^2. Note: The 4 lowest bits
// of a pure Y^2 would be zeros, but they are shifted out to preserve as much
// precision as possible. Within the 64 bits of YlowSqr, bits 18 to 31 will
// be used as a hash index to facilitate later searches.
//
// Also note: The trial edge "Y" may have more digits than the 15+ decimal
// digit precision in a "double" variable. The full precision is not needed to
// verify a solution. If a solution is found and the exact value of "Y" is too
// big for precision in a "double" variable, it would be relatively easy to
// modify the code to provide this higher precision.
//
// Finally, the Y's are split into those that are evenly divisible by 16 and
// those that are not. (Most will wind up in the 2nd group). To form a brick
// with integer dimensions, at least one side must come from the divide by 16
// group. Since no combination is possible if both sides are in the not div.
// by 16 group, this splitting will speed up later searching.
//
//*****************************************************************************
void AddEdge(double D1,double D2, double D3) {
// Note: "D2" will be > "D1", and all three
// parameters will be odd integers < 1.0E12
long long unsigned YlowSqr;
unsigned YsqrHash; // Bits 18 to 31 of YlowSqr
unsigned link, Yindex, Div16Test;
double p, q;
double Yedge;
long long unsigned Puns, Quns, Kuns; // Used to construct YlowSqr
// In terms of magnitude, "p", "q", and "k"
// will be < 1.0E12 as long as Xodd < 1.0E12
// Solve: P - Q = D1, P + Q = D2 for P and Q
p = (D1 + D2) / 2.0; // Add the above equations, and divide by 2.
q = p - D1; // Note that either p or q will be even
Yedge = 2.0 * p * q * D3; // Note: Good for about 15 decimal digits.
// If higher precision is ever needed, put the
// calcs here. (Might be able to get by with
// "long double" instead "double" for Yedge.)
// Next, get 64 low bits of the above.
Puns = p; // Convert p, q, and k to 64 bit integers
Quns = q;
if (Puns & 1) // Get rid of the trailing "0" bit that will
Quns >>= 1; // exist in either "P" or "Q"
else
Puns >>= 1;
Kuns = D3; // YlowSqr is just used for the look-up system
YlowSqr = Puns * Quns * Kuns; // Thus don't need to multiply by 2. Any
// overflow in YlowSqr can be ignored as we
// only need the lowest 64 bits.
YlowSqr *= YlowSqr; // Square it.
YsqrHash = YlowSqr; // Gets lowest 32 bits.
Div16Test = YsqrHash; // Will use to see if Yedge is divisible by 16
YsqrHash >>= 18; // Just use bits 18 to 31 for the hash index.
// Use hash and link list to see if
// "Y" is already in the table
for (link = HashHead[YsqrHash]; link; link = HashLink[link]) {
if (YlowSqr != YsqrLow[link]) // Most non matches will
continue; // occur here.
if (YforTrials[link] == Yedge) // Rarely needed.
break;
}
// Note: If "Yedge" is evenly divisible by 16, then it has at
// least 4 trailing binary zeroes. If Y/16 is an integer, then the
// expression
// "YlowSqr = Puns * Quns * Kuns;" (before squaring)
// will have at least 2 trailing zeroes. (2 were shifted/omited)
// If Y/16 is even, then after squaring YlowSqr will have at
// least 4 trailing zeroes. If we "AND" this with "15"
// (Binary 1111), then the result will be true if and only if 16
// does not divide Y evenly.
if (!link) { // If not found then add new edge to list
if (Div16Test & 15) Yindex = ++NotDiv16; // Loc. in Ylist is determined
else Yindex = --Div16; // by divisibility by 16.
YforTrials[Yindex] = Yedge; // Add y to candidate list
YsqrLow[Yindex] = YlowSqr; // Ditto for 64 low bits of Y^2.
HashLink[Yindex] = HashHead[YsqrHash];
HashHead[YsqrHash] = Yindex; // Insert in link list
} // End if not found
} // End AddEdge
//*****************************************************************************
//
// Search Edges
//
// This routine tries all valid combinations in the trial Y edges list to see
// if the sum of the squares of any two members matches any 3rd value in the
// list. If a match is found, then the integer brick problem has been solved.
//
// There may be many thousands of members in this list. All efforts have been
// used to speed up this search. The hash look-up system avoids the inner loop
// of what would otherwise look like:
//
// For each member of the group
// For each member of the group
// Add the squares of these two numbers together
// For each member of the group
// See if the square of this third member matches the above sum
// (If yes then output solution)
// Repeat third loop
// Repeat second loop
// Repeat first loop
//
//*****************************************************************************
void SearchEdges(void) {
unsigned i, j;
unsigned link, hash;
long long unsigned y1LowBits, SumLowBits;
double hypotenuse; // Used for external solutions where
// Xodd is less than 5000
for (i = Div16; i < 25000; i++) { // A brick requires at least 1 member
y1LowBits = YsqrLow[i]; // from this group.
for (j = i + 1; j <= NotDiv16; j++) {
// If in early stages of the program, then output external
// solutions. The 5,000 can be arbitrarily increased, but will
// rapidly run into precision problems. Note that precision
// calculations are only applied to potential complete solutions
// that include the internal diagonal.
if (Xodd < 5000.0) { // Compute hypotenuse of the two "Y"
hypotenuse = hypot(YforTrials[i], YforTrials[j]); // sides.
if (hypotenuse == round(hypotenuse)) { // Calcs involve the "even"
ExtCount++; // face of the brick.
printf("External sol. Nbr. %d at X = %'.0f Y = %'.0f Z = %'.0f\n",
ExtCount, Xodd, YforTrials[i], YforTrials[j]);
}
} // End of "If in early stages"
// If the sum of any two trial Yedge^2 entrees is equal to
// some other entry in the table, then a solution to the
// integer brick problem may have been found. Get this sum
// and then search the list to see if there is a match.
SumLowBits = YsqrLow[j] + y1LowBits;
link = SumLowBits; // Construct hashlink head for this
link >>= 18; // sum
for (link = HashHead[link]; link; link = HashLink[link]) {
if (SumLowBits != YsqrLow[link])
continue; // No solution.
// The following code is only used if a 64 bit match has been found.
// If the sum of the squares of any two sides is approximately the
// same as the square of some 3rd value, then a solution may have been
// found. If very high precision were available, then an exact match
// would be used. The "1.0E-12" limit might allow a "false positive"
// to be reported, but the relatively wide limit is used to prevent
// precision errors from skipping over a real solution.
if (fabs(hypot(YforTrials[i], YforTrials[j])/YforTrials[link] - 1.0)
< 1.0E-12) {
printf("\nPossible solution at X = %'.0f Y =%'.0f Z =%'.0f\n",
Xodd, YforTrials[i], YforTrials[j]);
printf("with the diagonal on the \"Even Side\" = %'.0f\n",
YforTrials[link]);
hash = 0;
while (hash != 3) { // Pause the program so the
puts("Enter \"3\" to continue"); // result can be observed.
gets(Databuff);
hash = atoi(Databuff);
}
} // End if low matches
} // End for link
} // End for j
} // End for i
// Reset HashHead array back to zeros.
for (i = Div16; i <= NotDiv16; i++) {
hash = YsqrLow[i];
hash >>= 18;
HashHead[hash] = 0;
}
} // End SearchEdges
Return to the main Integer Brick Page
Return to Durango Bill's Home Page
Web page generated via Sea Monkey's Composer HTML editor
within a Linux Cinnamon Mint 18 operating system.
(Goodbye Microsoft)