I extended the wheel factorization to a basis of 2, 3, 5, and 7 and added user input validation and output display in the following two routines. The text listings come from Takumako's CcEditor. See https://cclinker.web.fc2.com/
WFSUB
E+1→E↵
B→Z[E]↵
Do↵
D→A↵
Z[E+11]+1→Z[E+11]↵
A÷B→D↵
LpWhile Frac(D)=0↵
Int(√(A))→C↵
Return↵
FACTOR
0→DimZ↵
22→DimZ↵
"NUMBER"?→F↵
If F<1 Or F≧1x₁₀10↵
Then "NUMBER□MUST□BE□□≧1□ And <1x₁₀10"↵
Stop↵
IfEnd↵
If F≠Int(F)↵
Then "NUMBER□MUST□BE□□AN□INTEGER"↵
Stop↵
IfEnd↵
For 1→E To 22↵
0→Z[E]↵
Next↵
1→Z[1]↵
1→Z[12]↵
1→E↵
F→A↵
Int(√(A))→C↵
2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
3→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
5→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
7→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
11→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
While 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+8→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+8→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+6→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+4→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+10→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+2→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
B+10→B:A÷B→D:Frac(D)=0⇒Prog "WFSUB":B>C⇒Goto 1↵
WhileEnd↵
Lbl 1↵
If A>1↵
Then E+1→E↵
A→Z[E]↵
1→A↵
1→Z[E+11]↵
IfEnd↵
Int(E÷3)→D↵
E-3ｘD>0⇒D+1→D↵
1→C↵
Lbl 2↵
Cls↵
Locate 1,1,F↵
Locate 12,1,C↵
Locate 13,1,":"↵
Locate 14,1,D↵
3ｘ(C-1)+1→B↵
Locate 1,2,Z[B]↵
Locate 11,2,"^("↵
Locate 13,2,Z[B+11]↵
Locate 16,2,")"↵
If B+1≦E↵
Then Locate 1,3,Z[B+1]↵
Locate 11,3,"^("↵
Locate 13,3,Z[B+12]↵
Locate 16,3,")"↵
IfEnd↵
If B+2≦E↵
Then Locate 1,4,Z[B+2]↵
Locate 11,4,"^("↵
Locate 13,4,Z[B+13]↵
Locate 16,4,")"↵
IfEnd↵
While 1↵
Getkey→A↵
If A=34 Or A=73↵
Then Cls↵
"DONE"↵
Stop↵
IfEnd↵
If A=84 Or A=86 Or A=77 Or A=47↵
Then C+1→C↵
C>D⇒1→C↵
Goto 2↵
IfEnd↵
If A=83 Or A=85 Or A=67↵
Then C-1→C↵
C<1⇒D→C↵
Goto 2↵
IfEnd↵
WhileEnd↵
FACTOR can factorize 6666666667 in 20 seconds on the fx-5800P, but this number has some small factors which makes the trial division go more quickly. The worst-case input is 9999999769, which takes 18 minutes and 22 seconds to "factor"; this is the largest prime number less than 1x₁₀10.
I wrote the following test code in C++, and my computer is more than halfway done testing the factorization routine for all inputs on the range [1,9999999999].
#include <cinttypes>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <string>
#include <stdexcept>
/**
* Calculates the integer square root of a 64 bit unsigned number.
* From http://www.codecodex.com/wiki/Calculate_an_integer_square_root
* with stylistic modifications and extension to 64 bits.
*/
uint64_t isqrt64(const uint64_t &number)
{
uint64_t root = 0ull;
uint64_t remainder = number;
uint64_t place = 0x4000000000000000ull;
while (place > remainder) place >>= 2u;
while (place != 0ull)
{
if (remainder >= (root + place))
{
remainder -= root + place;
root += place << 1u;
}
root >>= 1u;
place >>= 2u;
}
return root;
}
/**
* Test function.
* Returns true if number is prime and false otherwise.
*/
bool isPrime(const uint64_t &number)
{
if (number == 0ull) return false;
if ((number == 1ull) | (number == 2ull)) return true;
if ((number & 1ull) == 0ull) return false;
uint64_t max_prime = isqrt64(number) + 1ull;
for (uint64_t test = 3ull; test <= max_prime; test += 2ull)
{
if (number % test == 0ull) return false;
}
return true;
}
/**
* Simulating the variables in the Casio fx-5800P with
* globals. To maintain similarity, Z_[0] is not used.
*/
double A;
double B;
double C;
double D;
double E;
double F;
double G;
double H;
double I;
double J;
double K;
double L;
double M;
double N;
double O;
double P;
double Q;
double R;
double S;
double T;
double U;
double V;
double W;
double X;
double Y;
double Z;
double Z_[23];
/**
* Implements the Casio fx-5800P Frac( function using
* the C++ standard library.
*/
inline double Frac(const double &number)
{
double scratch;
return std::modf(number, &scratch);
}
/**
* Subroutine of WHEELFAC()
* Call with B a known factor of A.
* Pulls all powers of B out of A and stores
* B as a new factor as specified below:
*
* 1. E holds number of factors (ignoring powers)
* and must lie on the range [1,11].
*
* a. E is limited to 11 because
* 1 * 2 * 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 * 29
* is the largest set of unique factors whose
* product is less than or equal to 10000000000, the
* largest number in the range represented as exact
* integers in the fx-5800P.
*
* 2. Z_[E] holds B.
*
* 3. Z_[E + 11] holds the highest power of B
* that is a factor of A.
*
* WFSUB modifies A, C, D, E, and Z_[1] through Z_[22].
* It uses but does not modify B.
*/
void WFSUB(void)
{
E = E + 1.0;
if (E > 11.0) { std::printf("\nERROR: ran out of factor storage space\n\n"); A = 1.0; return; }
Z_[static_cast<uint16_t>(E)] = B;
do
{
A = D;
Z_[static_cast<uint16_t>(E + 11.0)] = Z_[static_cast<uint16_t>(E + 11.0)] + 1.0;
D = A / B;
} while (Frac(D) == 0.0);
C = std::trunc(std::sqrt(A));
return;
}
/**
* Wheel factorization with a basis of 2, 3, 5, and 7.
* Stores factors of A in Z_[1] through Z_[22] per the
* comment above WFSUB.
* Modifies, A, B, C, D, E, and Z_[1] through Z_[22].
* Calls WFSUB.
*/
void WHEELFAC(void)
{
for (E = 1.0; E <= 22.0; E = E + 1.0) Z_[static_cast<uint16_t>(E)] = 0.0;
Z_[1] = 1.0;
Z_[12] = 1.0;
E = 1.0;
C = std::trunc(std::sqrt(A));
B = 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = 3.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = 5.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = 7.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = 11.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
while (true)
{
B = B + 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 4.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 4.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 6.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 6.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 4.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 4.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 6.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 6.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 6.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 4.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 6.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 4.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 6.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 8.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 4.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 4.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 4.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 8.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 6.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 4.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 6.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 4.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 6.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 6.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 6.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 4.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 4.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 6.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 6.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 4.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 4.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 10.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 2.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
B = B + 10.0; D = A / B; if (Frac(D) == 0.0) WFSUB(); if (B > C) goto one;
}
one:
if (A > 1.0)
{
E = E + 1.0;
Z_[static_cast<uint16_t>(E)] = A;
A = 1.0;
Z_[static_cast<uint16_t>(E + 11.0)] = 1.0;
}
return;
}
/**
* Test the WHEELFAC() routine.
* Take and validate user input to determine the range of numbers
* used during testing.
* The stride argument facilitates testing a large number range
* by running multiple instances of this program in parallel.
*/
int main(int argc, char *argv[])
{
std::printf("cmdline arguments: [start] [end] [stride]\n");
std::printf("all arguments required. must be positive integers.\n\n");
if (argc != 4)
{
std::printf("ERROR: must run with 3 arguments.\n\n");
return 0;
}
int64_t istart;
try { istart = std::stoll(argv[1]); }
catch (const std::out_of_range &e)
{
std::printf("ERROR: start argument exceeds range of int64\n");
std::printf(" pick start on range [1,9999999999]\n\n");
return 0;
}
catch (const std::invalid_argument &e)
{
std::printf("ERROR: could not convert start argument to integer\n\n");
return 0;
}
if ((istart < 1ll) | (istart > 9999999999ll))
{
std::printf("ERROR: start argument out of allowed range [1,9999999999]\n\n");
return 0;
}
int64_t iend;
try { iend = std::stoll(argv[2]); }
catch (const std::out_of_range &e)
{
std::printf("ERROR: end argument exceeds range of int64\n");
std::printf(" pick end on range [1,9999999999]\n\n");
return 0;
}
catch (const std::invalid_argument &e)
{
std::printf("ERROR: could not convert end argument to integer\n\n");
return 0;
}
if ((iend < 1ll) | (iend > 9999999999ll))
{
std::printf("ERROR: end argument out of allowed range [1,9999999999]\n\n");
return 0;
}
if (iend < istart)
{
std::printf("ERROR: end argument %" PRIi64 " less than start argument %" PRIi64 "\n\n", iend, istart);
return 0;
}
int64_t istride;
try { istride = std::stoll(argv[3]); }
catch (const std::out_of_range &e)
{
std::printf("ERROR: stride argument exceeds range of int64\n\n");
return 0;
}
catch (const std::invalid_argument &e)
{
std::printf("ERROR: could not convert stride argument to integer\n\n");
return 0;
}
if (istride < 1ll)
{
std::printf("ERROR: stride argument must be positive\n\n");
return 0;
}
if (istride > (iend - istart))
{
std::printf("WARNING: stride = %" PRIi64 " exceeds (end - start) = %" PRIi64 "\n", istride, iend - istart);
std::printf(" only start = %" PRIi64 " will be evaluated.\n\n", istart);
}
uint64_t total_evals = 1ull + (iend - istart) / istride;
uint64_t n_evals = 0ull;
uint64_t print_count = 0ull;
for (int64_t i = istart; i <= iend; i += istride)
{
double i_dbl = static_cast<double>(i);
A = i_dbl;
WHEELFAC();
double product = 1.0;
for (uint32_t j = 1u; j <= 11u; j += 1u)
{
if (Z_[j] == 0.0) break;
if (Z_[j] < 0.0) std::printf("\nERROR: negative factor %.0f for input %.0f\n\n", Z_[j], i_dbl);
if (Z_[j] != std::trunc(Z_[j])) std::printf("\nERROR: non-integer factor %.16f for input %.0f\n\n", Z_[j], i_dbl);
if (!isPrime(static_cast<uint64_t>(Z_[j]))) std::printf("\nERROR: composite factor %.0f for input %.0f\n\n", Z_[j], i_dbl);
for (double k = Z_[j + 11u]; k > 0.0; k = k - 1.0) product = product * Z_[j];
}
if (product != i_dbl)
{
std::printf("\nERROR: incorrect factorization for input %.0f\n\n", i_dbl);
for (uint32_t j = 1u; j <= 11u; j += 1u)
{
if (Z_[j] == 0.0) break;
std::printf("%.0f^%.0f\n", Z_[j], Z_[j + 11u]);
}
}
bool sorted = true;
for (uint32_t j = 1u; j <= 10u; j += 1u)
{
if (Z_[j] == 0.0) break;
if ((Z_[j + 1u] != 0.0) & (Z_[j + 1u] <= Z_[j]))
{
sorted = false;
break;
}
}
if (!sorted)
{
std::printf("\nERROR: factors not sorted for input %.0f\n\n", i_dbl);
for (uint32_t j = 1u; j <= 11u; j += 1u)
{
if (Z_[j] == 0.0) break;
std::printf("%.0f^%.0f\n", Z_[j], Z_[j + 11u]);
}
}
n_evals++;
print_count++;
if (print_count >= 10000000ull)
{
print_count = 0ull;
std::printf("progress: %" PRIu64 " / %" PRIu64 "\n", n_evals, total_evals);
}
}
std::printf("\nIf no errors above, it worked.\n");
std::printf("start = %" PRIi64 ", end = %" PRIi64 ", stride = %" PRIi64 "\n\n", istart, iend, istride);
return 0;
}
The code above, and program listings with notes and variable descriptions can be found at https://slugrustle.github.io/ and https://github.com/s.../fx-5800P_progs.
Edit 2018-12-12: My computer finished running the test code on 1, 2, ..., 9999999999. The algorithm worked correctly.
Edited by slugrustle, 12 December 2018 - 03:17 PM.
