Sequence of primorial primes
The sequence of primorial primes is given as the increasing values of n where primorial(n) ± 1 is prime.
Noting that the n'th primorial is defined as the multiplication of the smallest n primes, the sequence is of the number of primes, in order that when multiplied together is one-off being a prime number itself.
- Task
Generate and show here the first ten values of the sequence.
- Optional extended task
Show the first twenty members of the series.
- Notes
- This task asks for the primorial indices that create the final primorial prime numbers, so there should be no ten-or-more digit numbers in the program output (although extended precision integers will be needed for intermediate results).
- There is some confusion in the references, but for the purposes of this task the sequence begins with n = 1.
- Cf.
- Primorial prime Wikipedia.
- Primorial prime from The Prime Glossary.
- Sequence A088411 from The On-Line Encyclopedia of Integer Sequences
- Related tasks
C
<lang c>
- include <gmp.h>
int main(void) {
mpz_t p, s; mpz_init_set_ui(p, 1); mpz_init_set_ui(s, 1);
for (int n = 1, i = 0; i < 20; n++) {
mpz_nextprime(s, s);
mpz_mul(p, p, s);
mpz_add_ui(p, p, 1);
if (mpz_probab_prime_p(p, 25)) {
mpz_sub_ui(p, p, 1);
gmp_printf("%d\n", n);
i++;
continue;
}
mpz_sub_ui(p, p, 2);
if (mpz_probab_prime_p(p, 25)) {
mpz_add_ui(p, p, 1);
gmp_printf("%d\n", n);
i++;
continue;
}
mpz_add_ui(p, p, 1); }
mpz_clear(s); mpz_clear(p);
} </lang>
- Output:
1 2 3 4 5 6 11 13 24 66 68 75 167 171 172 287 310 352 384 457
Clojure
- Uses Java Interoop to use Java for 1) Prime Test and 2) Prime sequence generation
<lang lisp> (ns example
(:gen-class))
- Lazy Sequence of primes (starting with number 2)
(def primes (iterate #(.nextProbablePrime %) (biginteger 2)))
(defn primorial-prime? [v]
" Test if value is a primorial prime "
(let [a (biginteger (inc v))
b (biginteger (dec v))]
(or (.isProbablePrime a 16)
(.isProbablePrime b 16))))
- Generate indexes for first 20 primorial primes
(println (take 20 (keep-indexed ; take the first 20
#(if (primorial-prime? %2) (inc %1)) ; filters out non-primorials, passing on the index + 1 (since sequence begins with 1 (not 0)
(reductions *' primes)))) ; computes the lazy sequence of product of 1 prime, 2 primes, 3 primes, etc.
</lang>
- Output:
(1 2 3 4 5 6 11 13 24 66 68 75 167 171 172 287 310 352 384 457)
EchoLisp
<lang lisp> (lib 'timer) ;; for (every (proc t) interval) (lib 'bigint)
- memoize primorial
(define p1000 (cons 1 (primes 1000))) ; remember first 1000 primes (define (primorial n)
(if(zero? n) 1 (* (list-ref p1000 n) (primorial (1- n)))))
(remember 'primorial)
(define N 0)
- search one at a time
(define (search t) ;; time parameter, set by (every), not used (set! N (1+ N)) (if (or (prime? (1+ (primorial N))) (prime? (1- (primorial N)))) (writeln 'HIT N )
(writeln N (date->time-string (current-date )))))
</lang>
- Output:
<lang lisp>
- run in batch mode to make the browser happy
- search next N every 2000 msec
(every 2000 search)
→ HIT 1 HIT 2 HIT 3 HIT 4 HIT 5 HIT 6 7 "13:47:03" HIT 11 HIT 13 HIT 24 HIT 66 HIT 68 HIT 75 HIT 167 HIT 171 HIT 172 HIT 287 HIT 310 HIT 352 HIT 384 455 "14:07:08" 456 "14:07:12" HIT 457 458 "14:07:25" 459 "14:07:29"
</lang>
Fortran
The Plan
First comes a probe of the problem, using ordinary arithmetic and some routines written for other problems. Module PRIMEBAG can be found in Extensible_prime_generator#Fortran and the basic version of module BIGNUMBERS in Primorial_numbers#Fortran <lang Fortran> PROGRAM PRIMORIALP !Simple enough, with some assistants.
USE PRIMEBAG !Some prime numbers are wanted.
USE BIGNUMBERS !Just so.
TYPE(BIGNUM) B !I'll have one.
INTEGER MAXF !Largest factor to consider by direct division.
PARAMETER (MAXF = 18000000) !Some determination.
INTEGER I !Step stuff.
INTEGER FU,FD !Found factors.
INTEGER NHIT,HIT(666) !A little list.
CHARACTER*4 WOT !A remark.
CHARACTER*66 ALINE !A scratchpad.
REAL T0,T1 !In memory of lost time.
MSG = 6 !Standard output.
WRITE (MSG,1) BIGLIMIT,BIGBASE,HUGE(I) !Announce.
1 FORMAT ('Calculates primorial "primes"',/,
1 "A primorial prime is a value N such that",/,
2 " Primorial(N) - 1 is prime, OR",/,
3 " Primorial(N) + 1 is prime, or both.",/,
4 "and Primorial(N) is the product of the first N prime numbers.",/
5 "Working with up to ",I0," digits in base ",I0,"."/
6 "The integer limit is ",I0,/)
c CALL PREPARE PRIMES !First, catch your rabbit. Via ERATOSTHENES.
IF (.NOT.GRASPPRIMEBAG(66)) STOP "Gan't grab my file!" !Attempt in hope.
WRITE (MSG,2)
2 FORMAT ("Primorial#",3X,"Approx.",8X," -1 Factor +1 Factor Hit")
Commence prime mashing.
100 NHIT = 0 !My list is empty.
B.LAST = 1 !Begin at the beginning.
B.DIGIT(1) = 1 !With one. The same, whatever BIGBASE.
CALL CPU_TIME(T0) !Start the timing.
DO I = 1,30 !69 !Step along the primorials.
CALL BIGMULTN(B,PRIME(I)) !Multiply by the next prime.
c WRITE (MSG,101) I,PRIME(I),I,B.DIGIT(B.LAST:1:-1) !Digits in Arabic/Hindu order.
101 FORMAT ("Prime(",I0,") = ",I0,", Primorial(",I0,") = ", !For a possibly multi-BIGBASE sequence.
1 I0,9I<BIGORDER>.<BIGORDER>,/,(10I<BIGORDER>.<BIGORDER>)) !The first without leading zero digits.
FU = -1 !No factor for up one.
FD = -1 !No factor for down one.
CALL BIGADDN(B,+1) !Go up one.
FU = BIGFACTOR(B,MAXF) !Find a factor, maybe.
CALL BIGADDN(B,-2) !Now test down one.
IF (FU.NE.1) FD = BIGFACTOR(B,MAXF) !But only if FU didn't report "prime".
IF (FU.EQ.1 .OR. FD.EQ.1) THEN !Since if either candidate is a prime,
WOT = "Yes!" !Then a hit.
NHIT = NHIT + 1 !So count up a success.
HIT(NHIT) = I !And append to my list.
ELSE IF (FU.GT.1 .AND. FD.GT.1) THEN !But if both have factors,
WOT = "No." !Then definitely not a hit.
ELSE !Otherwise,
WOT = "?" !I can't decide.
END IF !So much for that candidate.
CALL BIGADDN(B,1) !Recover the original primorial value.
WRITE (ALINE,102) I,BIGVALUE(B),FD,FU,WOT !Prepare a report.
102 FORMAT (I10,1PE18.10,I10,I10,1X,A) !A table.
IF (FD.EQ.-1) ALINE(37:38) = "" !Wasn't looked for, so no remark.
IF (FD.EQ. 0) ALINE(38:38) = "?" !Recode a zero.
IF (FU.EQ. 0) ALINE(48:48) = "?" !Since it represents "Don't know".
WRITE (MSG,"(A)") ALINE !Show the report.
END DO !On to the next prime.
CALL CPU_TIME(T1) !Completed the run.
Cast forth some pearls.
WRITE (MSG,201) HIT(1:NHIT) !The list.
201 FORMAT (/,"Hit list: ",I0,666(",",I0:)) !Don't actually expect so many.
WRITE (MSG,*) "CPU time:",T1 - T0 !The cost.
END !So much for that.
</lang> Output, somewhat edited, because the limit on direct factorisation is soon exceeded: five lines such as "Passed the limit of 18000000 with 18000041, but Sqrt(n) = 0.43849291E+11:" and similar have been removed, corresponding to the five "?" reports.
Calculates primorial "primes"
A primorial prime is a value N such that
Primorial(N) - 1 is prime, OR
Primorial(N) + 1 is prime, or both.
and Primorial(N) is the product of the first N prime numbers.
Working with up to 8888 digits in base 10.
The integer limit is 2147483647
Primorial# Approx. -1 Factor +1 Factor Hit
1 2.0000000000E+00 1 Yes!
2 6.0000000000E+00 1 Yes!
3 3.0000000000E+01 1 Yes!
4 2.1000000000E+02 1 Yes!
5 2.3100000000E+03 1 Yes!
6 3.0030000000E+04 1 59 Yes!
7 5.1051000000E+05 61 19 No.
8 9.6996900000E+06 53 347 No.
9 2.2309287000E+08 37 317 No.
10 6.4696932300E+09 79 331 No.
11 2.0056049013E+11 1 Yes!
12 7.4207381348E+12 229 181 No.
13 3.0425026353E+14 1 61 Yes!
14 1.3082761332E+16 141269 167 No.
15 6.1488978259E+17 191 953 No.
16 3.2589158477E+19 87337 73 No.
17 1.9227603502E+21 ? 277 ?
18 1.1728838136E+23 1193 223 No.
19 7.8583215511E+24 163 ? ?
20 5.5794083013E+26 ? 1063 ?
21 4.0729680599E+28 313 2521 No.
22 3.2176447673E+30 163 22093 No.
23 2.6706451569E+32 139 265739 No.
24 2.3768741896E+34 ? 131 ?
25 2.3055679639E+36 66683 2336993 No.
26 2.3286236436E+38 ? 960703 ?
27 2.3984823529E+40 15649 2297 No.
28 2.5663761176E+42 17515703 149 No.
29 2.7973499682E+44 719 334507 No.
30 3.1610054640E+46 295201 5122427 No.
Hit list: 1,2,3,4,5,6,11,13
CPU time: 2.968750
The direct factorisation assault was extended up to 18,000,000 because one factor was 17,515,703 and it was pleasing to have a definite result. However there is not much hope for this approach because the 32-bit integer limit is exceeded after Primorial(9) and the 64-bit limit after Primorial(15), so definitive factorisation in a reasonable time is out of reach. Already, bignumber arithmetic has been introduced to calculate this far, and only the first eight candidates have been found when the demand is for at least ten.
Although there would be no difficulty in principle in a Fortran compiler accepting syntax such as INTEGER*1000, this sort of ability is not required in the modern standard. Fortran II of 1957 on a decimal computer such as the IBM1620 allowed a control card to state how many digits were to be used for fixed-point numbers (i.e. integers) but although the hardware delivered integer arithmetic of arbitrary precision (up to memory limits), thousand-digit arithmetic probably was out of bounds, just as INTEGER*8 or perhaps INTEGER*16 is the upper limit for modern compilers. Thus, there is no "built-in" facility for arithmetic involving larger numbers available via a standard specification. So, one must devise appropriate routines.
Functions or Subroutines?
Rather than always employing in-line arithmetic expressions of ever greater complexity, defining a function offers multiple advantages. This might be via an "arithmetic statement" function if it can be written as a single expression, but more generally, a function can be named and defined as a separate unit for compilation, and then used within arithmetic expressions. This is unexceptional, and easily understood. Whatever the parameters and workings of the function, its result is returned in the computer's accumulator, or top-of-stack, or a standard register. Thus, an expression such as y:=x; would produce code such as Load x; Store y; whereas with an expression such as y:=F(x); instead of the Load there is an invocation of the function which would perform its calculation then return, whereupon the code Store y; would be executed just as before. This is easily extended for usage within an expression such as y:=y + 7*F(x) - 3; and so forth.
Within the function's definition would be a statement such as F = expression followed by a return (perhaps by reaching the end of the function's code) but more complex functions introduce additional considerations. It is often the case that the final result will be developed via multiple statements and conditional branches. For instance, <lang Fortran> FUNCTION SUM(A,N)
DIMENSION A(*)
SUM = 0 !Unconditionally, clear SUM.
DO I = 1,N !There may be less than one element to sum.
SUM = A(I) + SUM
END DO !Some early compilers executed a DO-loop at least once.
END</lang>
Prior to F90, there were compilers for which this did not work. If an assignment to SUM were to be regarded as an assignment to the accumulator, as above, then when that is not followed by a RETURN but by other statements, the value that had been saved in the accumulator will be lost by the accumulator's usage in the execution of the subsequent statements. Indeed, some compilers would not allow any usage of SUM within an expression. Thus, a local variable, perhaps called ASUM, would be defined and used instead, and then, before exiting the function (by whatever path) one must remember to have an assignment SUM = ASUM
This difficulty no longer exists, and whether or not a local variable is used, it is possible that the compiler will generate code that avoids the unnecessary shuffling of data from one location to another. Or perhaps not. With the introduction of CHARACTER*n variables, the amount of data being copied about can now be much larger than an accumulator or a data "word" or two. Consider an expression such as TEXT = UPPERCASE(TEXT) It would be pleasant to think that the function's work area would be the storage belonging to the recipient of the result of the function so that no copying from place to place need delay the computer's speedy completion of its assigned task...
Or, one could define SUBROUTINE UPPERCASE(TEXT) that produces its result in-place, and deal patiently with the requirement to convert a simple expression into a suitable sequence of subroutine calls.
Long multiplication
The method taught in primary school is adequate, but variations are possible. For instance, instead of writing out the multi-line scheme of partial products then adding them in a second pass to produce the final product, the terms can be added as they are produced by working down a column rather than across a row - see the tableau in BIGMULT. This however requires an accumulator able to hold not just a two-digit number, but the sum of many such numbers given that the numbers have many digits. If calculations are proceeding with a base close to or equal to the word size of the computer's arithmetic, this becomes less convenient because then carries must be propagated to prevent overflow while working down a column rather than once after the column sum is complete. The expectation here is that the base will be some power of ten, and on a computer using binary arithmetic the word size will not be fully used.
More interesting is the opportunity of working left-to-right rather than the usual right-to-left, which is to say that the digits of the product will be produced from the highest order down, except for the propagation of carries leftwards to higher order digits. This in turn means that when multiplying A by B with the result to appear in A, the product can be placed in A as it is developed, rather than being copies from some work area, and it also allows BIGMULT(A,A) to calculate A² - though a dedicated routine BIGSQUARE(A) could employ its own trickery.
Long division
The straightforward way to calculate B mod M when a MOD operation is unavailable is too calculate B - B/M*M (hoping that an enthusiastic compiler won't convert .../M*M to a no-operation) where the division is integer division that truncates any remainder. Some languages such as pl/i retain the fractional part (producing a floating-point number), but division of integers in Fortran produces an integer result via truncation. Directly translating this formula requires a divide, a multiply, and a subtract, so one instead turns to consider variations on division.
Proffesor Knuth
In The Art of Comuter Programming, volume 2 Seminumerical Algorithims Professor Knuth remarks "Here the ordinary pencil-and-paper method involves a certain amount of guesswork and ingenuity on the part of the person doing the division; we must either eliminate this guesswork from the algorithm or develop some theory to explain it more carefully."
As will be remembered from primary school, at each step the key lies in the choice of a value for Q, the trial quotient for subtracting Q*M from the working number, with the digits of Q shifted left to align with the current high-order digits of B. One makes an estimate from considering the leading digits of B and M, and annoyingly, misguesses: the paper soon becomes decorated with side calculations of 3*Q, 7*Q, etc. to avoid repeated mistakes. In base ten there can be only ten possible values for Q at each cycle, but in larger bases any such list becomes ridiculously long and anyway, the chance of a recurrence of a value for Q falls.
Professor Knuth remarks that most computer arithmetic units allow for a double-word value to be divided by a single-word value with the hardware typically employing a register pair, and after the division one register holds the result and the other the remainder. Unfortunately, high-level languages typically provide no access to this facility and so code for bignumber arithmetic is littered by the likes of Digit:=DD mod Base; Carry:=DD div Base; - two divisions where one should suffice. If one allowed mod to be represented by \ and div by /, it would be nice to have syntax something like Digit,Carry:=DD/\Base;
The plan then is to guess a value for Q by inspecting the first two digits of B and the first only of M by calculating DD = B.Digit(B.Last)*BigBase + B.Digit(B.Last - 1) - this is a double-digit value, where B.Last fingers the highest-order digit of B and thus (B.Last - 1) the second digit. Then calculate Q = DD/M.Digit(M.Last) Immediate difficulties arise: in base ten this might be 31/2 and thus Q is a gross over-estimate. The estimate can be restrained: Q = MIN(DD/M.Digit(M.Last),BigBase - 1) but this is not the only problem.
How good is this estimate? Not very. The difficulty is that M is approximated by a single digit only, and when that digit is small, its relative accuracy is low. A lead digit of one would represent M = 1000... and also M = 1999... which involves almost a factor of two, thus the accuracy of DD/M.Digit(M.Last) varies by a factor of two likewise. The accuracy for higher digits is much better, to the degree that Professor Knuth recommends multiplying both B and M by a factor D such that the lead digit of the revised M is as large as possible: D = BigBase/(M.Digit(M.Last) + 1), but even so, Q may still be over-estimated, though now only by one and sometimes by two. Further tests are required to reduce the frequency of one-too-many, and preclude two-too-many.
Equipped with a good estimate for Q, the next step is to subtract Q*M from the high-order digits of B that align with the digits of M. Because now the full value of M is employed, the result can be negative: overshoot. If so, add back one M to produce the reduced value of B, ready for the next attempt with M shifted over to again align with the high-order digits of B.
Not followed...
Although the intention of BIGMOD(B,M) is for the value of B to hold the result, changing the value of M is likely to be unexpected, even if it is later changed back to what it had been. M might even be a constant, protected from change. Similarly, the rescale for B may require another digit and storage for this may not be available. Copies to working variables could be made, but anyway, although small compared to the effort of BIGMOD, the rescaling of big numbers is not so trivial. Instead, recognising that in high-level languages there is no access to to the double/single features of the hardware mentioned above, one can rely on the availability of equal-sized division - that is, one word divided by one word, or double by double. Thus, there is the opportunity for DD to be the first two digits of B and MM the first two of M, and then Q = DD/MM The minimum accuracy of this calculation is always better than one digit since (in base ten) MM = 10 represents 100... to 109... and so no tricky adjustment checks are required. There is still the overflow issue with the required add back when using the full precision of M, but the overflow is always only by one extra M value and so the shift and add ploy discussed below is safe.
The intention is to calculate B mod M by subtracting from B various multiples of M until the result no longer exceeds M; this is the required remainder. This will be done by following the method of long division except without bothering to record the quotient, although in other use producing both the quotient and the remainder in one go may be helpful. To this end, BIGMOD(B,M) first checks for a single-digit M and for that case, invokes BIGMODN(B,M.Digit(1)): any further proceedings can be assured of at least two digits in M so that MM may be produced. Iteration proceeds to reduce B while B > M. Subroutine BIGNORM ensures that any leading zero digits produced (possibly followed by existing zero digits) will be trimmed off ready for BIGSIGN to compare the new B against M.
Consequences
Alas, DD/MM sometimes yields zero, as in 1001/11 (where DD/M.Digit(M.Last) would yield 10, reduced by the MIN to 9) or 9800/99, or 4100/588, and many other cases. One could pre-emptively subtract M once and then engage in the overflow recovery procedure, or, shift M right one digit place and reconsider. This means calculating Q = (DD*BIGBASE + B.DIGIT(L - 1))/MM and proceeding with that, again at a risk of going one too far when using M instead of its approximation MM. Even if the first two digits of M exceed the first two digits of B, this shift is always possible because B > M is assured via the earlier test. However, this recalculation involves a three-digit value (DD is already a two-digit value) and until now the working assumption has been for arithmetic of up to two-digits only. This means the size of the big base is limited. The difficulty is that the proper value for Q is anything from zero to BIGBASE - 1, and with division by a two-digit number MM, for the higher values of Q to emerge, DD must have a three-digit number. This means that calculating a provisional Q value (from DD divided by the lead digit of M) then considering Q*MM - (top three digits of B) to determine an improved estimate will still have trouble. A multi-digit divide isn't available (as via recursion) since this is the multi-digit divide scheme and it is in trouble. Still further possibilities open up by considering the reciprocal of MM (calculated once at the start) so that DD/MM can be calculated using multiplication, but enough: accept a severe limitation on the possible base of arithmetic. For various powers of ten, when using two's complement signed integers,
Word size: 16-bit 32-bit 64-bit Limit 3E+4 2E+9 9E+18 Base 10 = E1: b**4 b**9 b**18 Base 100 = E2: b**2 b**4 b**9 Base 1000 = E3: b b**3 b**6 Base 10000 = E4: b b**2 b**4 Base 100000 = E5: - b b**3 Base 1000000 = E6: - b b**3 Base 10000000 = E7: - b b**2 Base 100000000 = E8: - b b**2 Base 1000000000 = E9: - b b**2
Thus, when 64-bit signed integer arithmetic is available, the highest possible decimal base is one million if triple-digit numbers are contemplated, but if only double-digit numbers are needed, base one thousand million can be used, which is a reasonable fit in a 32-bit variable.
Since these proceedings are conducted in a high-level language, not machine code, floating-point arithmetic is available and could be used for the DD/MM estimates, with attention to rounding possible. However, integer arithmetic is exact while floating-point is not. If the digits were 32-bit integers, DD and MM would require 64-bit precision, and a 64-bit floating-point number does not offer that. Thus, there would still be constraints on the maximum utilisation of a digit's integer variable. So cheating won't help.
Overshoot
Anyone who has cranked the handle of a mechanical calculator such as the Odhner will recall a better procedure for when the accumulator goes negative and a loud "Ding!" sounds. Rather than adding back M to the accumulator then shifting one and proceeding at the next digit position, instead shift one then add rather than subtract M with each crank of the handle (hearing the crunching of digit cogs) until the accumulator again changes sign ("Ding!") to positive. Hearing these "dings" soon elicits automatic reversal of rotation. This is the equivalent of multiplying by nine through shifting one left to add once (i.e. ten times), then shifting back right to subtract once: two crank twirls, rather than adding nine times with nine twirls and no shifts. Naturally, Q*M is subtracted once (rather than subtracting M with each of Q turns of the crank handle), and this is done without multiplying out Q*M first.
This recovery from overshoot relies on having overflowed B by no more than one too many multiples of M, which has been ensured by the calculation of Q, Perhaps (in base ten) Q should have been 2.6 rather than 3; if so, shifting M along one place means that instead of 30*M, the equivalent of 26*M should have been used, and this can be attained by adding back 4*M. By this means the recovery also reduces B by at least one digit place. Reversing the overflow by adding back one M costs a multi-digit add without such progress.
Just as the subtraction of Q*M could make zero multiple leading digits of B, not just one, so also can an overflow lead to multiple leading zero digits after the recovery. The extreme case is for a long M number, with the overflow being by just one in the units position of M as currently aligned with B. The overflow will manifest as complemented numbers in the high-order digits, a sequence such as 99999123 (in base ten) with an outstanding borrow, or C = -1. Rather than shifting M over one place for the add back, a multi-place shift will enable a single add back to clear many leading digits, though the shifts must go no further than aligning the units digit of M and B since only integers are in consideration.
In 16-Bit Modern Microcomputers, G.W. Corsline remarks "... multiplication and division are somewhat difficult and involve rather esoteric solutions in the general case." (p193) and "... although straightforward, are not short, simple, or easy to understand." (p124).
Storage
Run-time allocation
As usual, the question "How long is a piece of string?" arises with regard to the arrays of digits to be manipulated. In old-style Fortran, one would define arrays with a size that seemed "big enough" and have an associated variable that counted the number of digits in use. Having two variables to be used in tandem was tiresome and invited mistakes: a possible ploy would be to reserve the first element of the array for the count, with the digits following. The EQUIVALENCE statement could be used to help in some ways (though this is not allowed for parameters to routines), somewhat as follows: <lang Fortran> INTEGER XLAST,XDIGIT(66),XFIELD(67)
EQUIVALENCE (XFIELD(1),XLAST) !The count at the start.
EQUIVALENCE (XFIELD(3),XDIGIT(1)) !Followed by XDIGIT(0)</lang>
This introduces three names, if in a systematic way... But such array indexing invites its own mistakes, especially when relying on lax bound checking. For example, an accidental assignment to XDIGIT(-1) would change XLAST.
With F90 came the ability to declare arrays with a lower bound other than one, and also to define compound data aggregates. Rather dubious code can be replaced by something like <lang Fortran> TYPE(BIGNUM) !Define an aggregate.
INTEGER LAST !Count of digits in use.
INTEGER DIGIT(many) !The first digit is number one.
END TYPE BIGNUM
TYPE(BIGNUM) X !I'll have one.</lang>
And in this case the X.LAST variable is not a part of the digit array and so can't be damaged via misindexing X.DIGIT as with the earlier attempt using EQUIVALENCE - providing that array bound checking is engaged. Thus, one gains the ability to use X in a way that is somewhat more self-documenting, but there is still the annoyance of deciding on how many digits will be enough, and further, the BIGNUM type allows for only a fixed size. A different size would entail a different type and thereby cause parameter mismatching - should that be checked...
A further feature of F90 enables the determination of storage needs at run time, and via the ALLOCATE statement, a variable may be defined that is indeed "big enough" for the expected need. First, change the specification of DIGIT in BIGNUM to something like INTEGER, ALLOCATABLE:: DIGIT(:) then after calculating a suitable value for MANY, put ALLOCATE(X.DIGIT(MANY)) By this means, different-sized variables all have the same type, and can be passed as parameters, etc. But before acclaiming this as another demonstration of the merit of modern features, consider the following code snippet from BIGMULTN, as produced by the Compaq F90/95 compiler. First, with fixed-size DIGIT arrays:
153: DO I = 1,B.LAST !Step through the digits, upwards powers. 004018E1 mov esi,dword ptr [ebx] 004018E3 mov dword ptr [ebp-8],esi 004018E6 mov edi,1 004018EB mov dword ptr [I (0046b9bc)],edi 004018F1 cmp esi,0 004018F4 jle BIGNUMBERS_mp_BIGMULTN+260h (00401adf) 154: D = B.DIGIT(I)*V + C !Grab a digit and apply the multiply. 004018FA cmp dword ptr [I (0046b9bc)],1 00401901 jl BIGNUMBERS_mp_BIGMULTN+90h (0040190f) 00401903 cmp dword ptr [I (0046b9bc)],8AEh 0040190D jle BIGNUMBERS_mp_BIGMULTN+99h (00401918) 0040190F xor eax,eax 00401911 mov dword ptr [ebp-34h],eax 00401914 dec eax 00401915 bound eax,qword ptr [ebp-34h] 00401918 imul edi,dword ptr [I (0046b9bc)],4 0040191F mov esi,dword ptr [V (0046b9b8)] 00401925 imul esi,dword ptr [ebx+edi] 00401929 add esi,dword ptr [C (0046b9c4)] 0040192F mov dword ptr [D (0046b9c0)],esi 155: B.DIGIT(I) = MOD(D,BIGBASE) !Place the resulting digit. 00401935 cmp dword ptr [I (0046b9bc)],1 0040193C jl BIGNUMBERS_mp_BIGMULTN+0CBh (0040194a) 0040193E cmp dword ptr [I (0046b9bc)],8AEh 00401948 jle BIGNUMBERS_mp_BIGMULTN+0D4h (00401953) 0040194A xor eax,eax 0040194C mov dword ptr [ebp-34h],eax 0040194F dec eax 00401950 bound eax,qword ptr [ebp-34h] 00401953 imul edi,dword ptr [I (0046b9bc)],4 0040195A mov esi,2710h 0040195F mov eax,dword ptr [D (0046b9c0)] 00401965 cdq 00401966 idiv eax,esi 00401968 mov eax,edx 0040196A mov dword ptr [ebx+edi],eax 156: C = D/BIGBASE !Agony! TWO divisions per step!!
With the ALLOCATE usage, the same source code leads to...
154: DO I = 1,B.LAST !Step through the digits, upwards powers. 00401A55 mov ebx,dword ptr [B] 00401A58 mov edx,dword ptr [ebx] 00401A5A mov dword ptr [ebp-8],edx 00401A5D mov eax,1 00401A62 mov dword ptr [I (0046cacc)],eax 00401A68 cmp edx,0 00401A6B jle BIGNUMBERS_mp_BIGMULTN+2E3h (00401ca2) 155: D = B.DIGIT(I)*V + C !Grab a digit and apply the multiply. 00401A71 mov ebx,dword ptr [B] 00401A74 mov esi,ebx 00401A76 mov edi,dword ptr [esi+20h] 00401A79 add edi,dword ptr [esi+18h] 00401A7C sub edi,1 00401A7F mov eax,dword ptr [I (0046cacc)] 00401A85 cmp eax,dword ptr [esi+20h] 00401A88 jl BIGNUMBERS_mp_BIGMULTN+0CFh (00401a8e) 00401A8A cmp eax,edi 00401A8C jle BIGNUMBERS_mp_BIGMULTN+0D8h (00401a97) 00401A8E xor edi,edi 00401A90 mov dword ptr [ebp-34h],edi 00401A93 dec edi 00401A94 bound edi,qword ptr [ebp-34h] 00401A97 imul eax,eax,4 00401A9A imul edx,dword ptr [esi+20h],4 00401A9E mov edi,dword ptr [esi+4] 00401AA1 sub edi,edx 00401AA3 mov edi,dword ptr [edi+eax] 00401AA6 imul edi,dword ptr [V (0046cac8)] 00401AAD add edi,dword ptr [C (0046cad4)] 00401AB3 lea esi,[esi+4] 00401AB6 mov dword ptr [D (0046cad0)],edi 156: B.DIGIT(I) = MOD(D,BIGBASE) !Place the resulting digit. 00401ABC mov ebx,dword ptr [B] 00401ABF mov ecx,ebx 00401AC1 mov eax,dword ptr [ecx+20h] 00401AC4 add eax,dword ptr [ecx+18h] 00401AC7 sub eax,1 00401ACC mov edx,dword ptr [I (0046cacc)] 00401AD2 cmp edx,dword ptr [ecx+20h] 00401AD5 jl BIGNUMBERS_mp_BIGMULTN+11Ch (00401adb) 00401AD7 cmp edx,eax 00401AD9 jle BIGNUMBERS_mp_BIGMULTN+125h (00401ae4) 00401ADB xor eax,eax 00401ADD mov dword ptr [ebp-34h],eax 00401AE0 dec eax 00401AE1 bound eax,qword ptr [ebp-34h] 00401AE4 imul esi,edx,4 00401AE7 imul edi,dword ptr [ecx+20h],4 00401AEB mov eax,dword ptr [ecx+4] 00401AEE sub eax,edi 00401AF0 mov edi,eax 00401AF2 mov ebx,2710h 00401AF7 mov eax,dword ptr [D (0046cad0)] 00401AFD cdq 00401AFE idiv eax,ebx 00401B00 mov eax,edx 00401B02 lea ecx,[ecx+4] 00401B05 mov dword ptr [edi+esi],eax 157: C = D/BIGBASE !Agony! TWO divisions per step!!
Evidently, every reference to an allocatable-size array involves quite a lot more code. A test run up to primorial 80 with run-time allocations of BIGLIMIT took ~4·42 seconds while fixed-size runs took ~4·13, a 7% increase. This is discouraging.
Simple array
In Fortran there has never been a problem with passing arrays of a different size to a subroutine (unlike in Pascal), so, representing the big number in a simple array (with DIGIT.LAST as the first element) would allow diferent sizes of a big number. One will just have to take care with the indexing and offsets, especially if indexing always starts with one. It so happens that type BIGNUM declared array DIGIT starting with DIGIT(1), so, taking advantage of F90's ability to allow indexing to start with zero enables B.LAST to be replaced by B(0) and B.DIGIT(i) to be replaced by B(i), though if the array indexing had started with zero, B(-1) could be reserved for the count. The type BIGNUM vanishes in favour of a simple INTEGER B(0:BIGLIMIT) - all this at a loss of self-documentation, though no fearsome EQUIVALENCE statements. The resulting test run takes about the same time at ~4·12 seconds, but, when the parameter declarations are changed from INTEGER B(0:BIGLIMIT) to INTEGER B(0:) to signify that the upper bound is flexible, the test run now takes ~4·5 seconds. Here is the same code from BIGMULTN for when the upper bound is specified:
153: DO I = 1,B(0) !Step through the digits, upwards powers. 004018CE mov esi,dword ptr [ebx] 004018D0 mov dword ptr [ebp-8],esi 004018D3 mov edi,1 004018D8 mov dword ptr [I (0046b9bc)],edi 004018DE cmp esi,0 004018E1 jle BIGNUMBERS_mp_BIGMULTN+260h (00401acc) 154: D = B(I)*V + C !Grab a digit and apply the multiply. 004018E7 cmp dword ptr [I (0046b9bc)],0 004018EE jl BIGNUMBERS_mp_BIGMULTN+90h (004018fc) 004018F0 cmp dword ptr [I (0046b9bc)],8AEh 004018FA jle BIGNUMBERS_mp_BIGMULTN+99h (00401905) 004018FC xor eax,eax 004018FE mov dword ptr [ebp-34h],eax 00401901 dec eax 00401902 bound eax,qword ptr [ebp-34h] 00401905 imul edi,dword ptr [I (0046b9bc)],4 0040190C mov esi,dword ptr [V (0046b9b8)] 00401912 imul esi,dword ptr [ebx+edi] 00401916 add esi,dword ptr [C (0046b9c4)] 0040191C mov dword ptr [D (0046b9c0)],esi 155: B(I) = MOD(D,BIGBASE) !Place the resulting digit. 00401922 cmp dword ptr [I (0046b9bc)],0 00401929 jl BIGNUMBERS_mp_BIGMULTN+0CBh (00401937) 0040192B cmp dword ptr [I (0046b9bc)],8AEh 00401935 jle BIGNUMBERS_mp_BIGMULTN+0D4h (00401940) 00401937 xor eax,eax 00401939 mov dword ptr [ebp-34h],eax 0040193C dec eax 0040193D bound eax,qword ptr [ebp-34h] 00401940 imul edi,dword ptr [I (0046b9bc)],4 00401947 mov esi,2710h 0040194C mov eax,dword ptr [D (0046b9c0)] 00401952 cdq 00401953 idiv eax,esi 00401955 mov eax,edx 00401957 mov dword ptr [ebx+edi],eax 156: C = D/BIGBASE !Agony! TWO divisions per step!!
It is much the same as for the non-allocated BIGNUM code. But when the specification is INTEGER B(0:) it becomes:
153: DO I = 1,B(0) !Step through the digits, upwards powers. 00402021 mov edi,dword ptr [ebx+1Ch] 00402024 add edi,dword ptr [ebx+14h] 00402027 sub edi,dword ptr [ebx+1Ch] 0040202A sub edi,1 0040202D mov eax,dword ptr [ebx] 0040202F xor esi,esi 00402031 cmp esi,edi 00402033 jle BIGNUMBERS_mp_BIGMULTN+0E1h (0040203e) 00402035 xor edi,edi 00402037 mov dword ptr [ebp-38h],edi 0040203A dec edi 0040203B bound edi,qword ptr [ebp-38h] 0040203E imul esi,dword ptr [ebx+18h] 00402042 mov edi,dword ptr [eax+esi] 00402045 mov dword ptr [ebp-8],edi 00402048 mov esi,1 0040204D mov dword ptr [I (0046d9bc)],esi 00402053 cmp edi,0 00402056 jle BIGNUMBERS_mp_BIGMULTN+327h (00402284) 154: D = B(I)*V + C !Grab a digit and apply the multiply. 0040205C mov eax,dword ptr [ebx+1Ch] 0040205F add eax,dword ptr [ebx+14h] 00402062 sub eax,dword ptr [ebx+1Ch] 00402065 sub eax,1 0040206A mov esi,dword ptr [ebx] 0040206C mov edi,dword ptr [I (0046d9bc)] 00402072 cmp edi,0 00402075 jl BIGNUMBERS_mp_BIGMULTN+11Eh (0040207b) 00402077 cmp edi,eax 00402079 jle BIGNUMBERS_mp_BIGMULTN+127h (00402084) 0040207B xor eax,eax 0040207D mov dword ptr [ebp-38h],eax 00402080 dec eax 00402081 bound eax,qword ptr [ebp-38h] 00402084 imul edi,dword ptr [ebx+18h] 00402088 mov esi,dword ptr [esi+edi] 0040208B imul esi,dword ptr [V (0046d9b8)] 00402092 add esi,dword ptr [C (0046d9c4)] 00402098 mov dword ptr [D (0046d9c0)],esi 155: B(I) = MOD(D,BIGBASE) !Place the resulting digit. 0040209E mov edi,dword ptr [ebx+1Ch] 004020A1 add edi,dword ptr [ebx+14h] 004020A4 sub edi,dword ptr [ebx+1Ch] 004020A7 sub edi,1 004020AA mov esi,dword ptr [ebx] 004020AC mov eax,dword ptr [I (0046d9bc)] 004020B2 cmp eax,0 004020B7 jl BIGNUMBERS_mp_BIGMULTN+160h (004020bd) 004020B9 cmp eax,edi 004020BB jle BIGNUMBERS_mp_BIGMULTN+169h (004020c6) 004020BD xor edi,edi 004020BF mov dword ptr [ebp-38h],edi 004020C2 dec edi 004020C3 bound edi,qword ptr [ebp-38h] 004020C6 imul eax,dword ptr [ebx+18h] 004020CA mov ecx,eax 004020CC mov edi,2710h 004020D1 mov eax,dword ptr [D (0046d9c0)] 004020D7 cdq 004020D8 idiv eax,edi 004020DA mov eax,edx 004020DC mov dword ptr [esi+ecx],eax 156: C = D/BIGBASE !Agony! TWO divisions per step!!
Accessing the secret additional parameter that gives the array parameter's actual upper bound for array bound checking takes more time. All this is a long way away from notions such as D = B.DIGIT(I)*V + C producing code like
Load I
Index B.DIGIT
Mult V
Add C
Store D
Fixed-size data aggregates
Since the BIGNUM type with fixed storage runs no slower but offers better documentation, the possibility of big numbers sized to fit is abandoned for demonstration purposes. Storage is plentiful these days, and can be expended with far less angst.
Source
Organisation
This uses the MODULE facilities of F90 to save on communication (though it seems that fixed-size work areas in COMMON may enable faster code) and although the DIGIT array could start with index zero so that the index value would correspond to the powers of the base, it starts at one so that the first digit, the units digit, is at index one. As ever, care is required with the various counting on digits.
Most of the source code is devoted to providing the BIGNUM facilities, since there is no built-in facility for this nor a universally-used library such as the "stdio.h" collection in C and similar languages. Via a great deal of additional syntax it is possible to extend F90 so as to provide the appearance of multi-precision arithmetic via the ordinary usage of expressions, as is exemplified for rational arithmetic in Arithmetic/Rational#Fortran. However, the "formula translation" that this invites will not work well for the likes of MOD(A**B,M) when the variables are hundred-digit numbers because A**B will overflow any possible storage scheme. The modulo arithmetic feature must be pushed into the calculation of A**B, and even if one has prepared a suitable routine such as BIGMODEXP, the language extension facilities are not going to recognise its suitability for calculating MOD(A**B,M) in a general context because they manifest "formula translation". One must employ explicit coding, as in BIGMRPRIME.
Support
<lang Fortran> MODULE BIGNUMBERS !Limited services: decimal integers, no negative numbers.
INTEGER BIGORDER !A limited attempt at generality.
PARAMETER (BIGORDER = 4) !This is the order of the base of the big number arithmetic.
INTEGER BIGBASE,BIGLIMIT !Sized thusly.
PARAMETER (BIGBASE = 10**BIGORDER, BIGLIMIT = 8888/BIGORDER) !Enough?
TYPE BIGNUM !So, a big number is simple.
INTEGER LAST !This many digits (of size BIGBASE) are in use.
INTEGER DIGIT(BIGLIMIT) !The digits, in ascending power order.
END TYPE BIGNUM !So much for that.
INTEGER BIGLEADN !Additional stuff for its rounding.
PARAMETER (BIGLEADN = 8/BIGORDER) !Sufficient digits to show. With a struggle.
INTEGER BIGLEAD(BIGLEADN + 1) !Followed by an exponent.
Collect some statistics on the working of BIGMRPRIME
INTEGER BIGMRTRIALS !How many trials, at most.
PARAMETER (BIGMRTRIALS = 6) !This might do.
INTEGER BIGMRCOUNT(BIGMRTRIALS) !BIGMRPRIME may not use all its trials.
DATA BIGMRCOUNT/BIGMRTRIALS*0/ !None so far.
CONTAINS !Now for some assistants.
SUBROUTINE BIGWRITE(F,B) !Show B.
INTEGER F !I/O unit number.
TYPE(BIGNUM) B !The number.
WRITE (F,1,ADVANCE="NO") B.DIGIT(B.LAST:1:-1) !Roll the digits in base ten.
1 FORMAT (665I<BIGORDER + 1>.<BIGORDER>) !Leading zeroes after the first digit.
END SUBROUTINE BIGWRITE !Simple, but messy.
INTEGER FUNCTION BIGSIGN(A,B) !Sign(A - B) returns -ve, 0, +ve. Not just -1, 0, +1.
TYPE(BIGNUM) A,B !The two numbers.
INTEGER L !A finger.
BIGSIGN = A.LAST - B.LAST !Compare the number of digits.
IF (BIGSIGN.EQ.0) THEN !The same?
DO L = A.LAST,1,-1 !Alas. Compare the digits themselves.
BIGSIGN = A.DIGIT(L) - B.DIGIT(L) !From the high-order, down.
IF (BIGSIGN.NE.0) RETURN !A difference yet?
END DO !Descend to the next pair of digits.
END IF !So much for shortcuts.
END FUNCTION BIGSIGN !Hopefully, not many digits will have to be compared.
SUBROUTINE BIGLOADN(B,N)!B:=N !Convert a normal number.
TYPE(BIGNUM) B !The multi-BIGBASE number.
INTEGER N !An ordinary number. Negatives need not apply.
INTEGER C !The carry.
B.LAST = 1 !Start with one digit.
C = N !Perhaps more than one will be needed.
1 B.DIGIT(B.LAST) = MOD(C,BIGBASE) !Place the digit.
C = C/BIGBASE !Lose that digit.
IF (C.GT.0) THEN !More remains?
B.LAST = B.LAST + 1 !Yes. Another digit is needed.
GO TO 1 !And try again.
END IF !Just in case N >= BIGBASE.
END SUBROUTINE BIGLOADN !Worry over negative numbers sometime later.
DOUBLE PRECISION FUNCTION BIGVALUE(B) !Convert back to an ordinary number. Limited range!
TYPE(BIGNUM) B !The mysterious big number.
REAL*8 F !The mantissa to come. "Infinity" starts at ~10**306.
INTEGER D !Counts BIGBASE digits.
INTEGER I !A stepper.
F = 0 !I'm not messing with negative numbers yet.
D = 0 !No digits assimilated.
DO I = B.LAST,MAX(1,B.LAST - 18/BIGORDER),-1 !Work from the high order down.
D = D + 1 !Another one in.
F = F*BIGBASE + B.DIGIT(I) !Thus.
END DO !Perhaps more to come.
BIGVALUE = F*DFLOAT(BIGBASE)**(B.LAST - D) !Beware floating-point overflow!
END FUNCTION BIGVALUE !An alternative would be BIGLOG(B).
SUBROUTINE BIGTASTE(B) !Represent B in BIGLEAD as a floating-point number.
Copies the first few digits to BIGLEAD, with rounding, followed by the base ten exponent.
TYPE(BIGNUM) B !The number to taste.
INTEGER L,IT !Fingers.
INTEGER D,E !A digit and an exponent.
BIGLEAD = 0 !Scrub, on general principles.
L = MIN(BIGLEADN,B.LAST) !I'm looking to taste the leading digits; too few?.
BIGLEAD(1:L) = B.DIGIT(B.LAST:B.LAST - L + 1:-1) !Reverse, to have normal order.
E = (B.LAST - 1)*BIGORDER !Convert from 10**BIGORDER to base 10.
D = B.DIGIT(B.LAST) !Grab the high-order digit.
DO WHILE(D.GT.0) !It is not zero..
E = E + 1 !So it is at least one base ten digit.
D = D/10 !Snip.
END DO !And perhaps there will be more.
D = 0 !We should consider rounding up.
IF (B.LAST.GT.BIGLEADN) THEN !Are there even more digits?
IT = L !Yes. This is now the low-order digit tasted.
IF (B.DIGIT(B.LAST - L).GE.BIGBASE/2) D = 1 !If the next digit is big enough.
C:DO WHILE (D.GT.0) !Spread the carry.
D = 0 !This one is used up.
BIGLEAD(IT) = BIGLEAD(IT) + 1 !Thusly.
IF (BIGLEAD(IT).GE.BIGBASE) THEN !But, maybe, overflow!
BIGLEAD(IT) = BIGLEAD(IT) - BIGBASE !Yes!
D = 1 !Reassert a carry.
IT = IT - 1 !Step back to the next recipient up.
IF (IT.LE.0) EXIT C !If there isn't one, quit!
END IF !So much for that overflow.
END DO C !Mostly, a carry doesn't propagate far.
END IF !So, no test for IT > 0 in a compound "while".
IF (D.NE.0) THEN !If a carry remains, the rounding propagated all the way up!
E = E + 1 !So, count another power of ten up.
BIGLEAD(1) = 1 !And thus maintain the 0.etc. style.
c BIGLEAD(2:) = 0 !These are already zero.
END IF !Enough carrying.
BIGLEAD(BIGLEADN + 1) = E !Place the exponent.
END SUBROUTINE BIGTASTE !That was messy.
SUBROUTINE BIGNORM(B) !Normalise B to avoide a horde of leading zero digits.
TYPE(BIGNUM) B !The number.
Can't rely on DO WHILE(B.LAST.GT.1 .AND. B.DIGIT(B.LAST)).NE.0) beause *both* parts *might* be evaluated. Bah.
10 IF (B.LAST.GT.1) THEN !If B has a single digit only, it can be zero.
IF (B.DIGIT(B.LAST).NE.0) RETURN !If it is not zero, we're done.
B.LAST = B.LAST - 1 !Otherwise, step down one digit,
GO TO 10 !And check afresh.
END IF !So much for normalisation. Rather like floating-point numbers.
END SUBROUTINE BIGNORM !In memory of Norman Kirk, Labour party leader and prime minister.
SUBROUTINE BIGADDN(B,N) !B:=B + N Add a (small) ordinary number to a big number.
TYPE(BIGNUM) B !The big number to be augmented.
INTEGER N !The addend. Should be smaller than the integer limit less BIGBASE.
INTEGER I !A stepper.
INTEGER C,D !Assistants for the arithmetic.
C = N !Start as if already in progress.
DO I = 1,B.LAST !Spread the carry to higher digits as needed.
D = B.DIGIT(I) + C !Thus.
IF (D.GE.0) THEN !Negative N might require a borrow.
B.DIGIT(I) = MOD(D,BIGBASE) !Not this time. Correct the digit.
C = D/BIGBASE !And calculate the carry to continue.
ELSE !Otherwise, borrow one from the next digit up.
D = D + BIGBASE !Thus.
B.DIGIT(I) = MOD(D,BIGBASE) !MOD for negative number can behave unexpectedly.
C = D/BIGBASE - 1 !Repay the borrow.
END IF !So much for DIGIT(I).
IF (C.EQ.0) RETURN !Finished already?
END DO !On to the next digit up.
Completed work with the original number of digits in the big number.
DO WHILE(C .NE. 0) !Now spread the last carry to further digits.
B.LAST = B.LAST + 1 !Up one more.
IF (B.LAST .GT. BIGLIMIT) STOP "Overflow by addition!" !Perhaps not.
B.DIGIT(B.LAST) = MOD(C,BIGBASE) !The digit.
C = C/BIGBASE !The carry may be large, if N is large.
END DO !So slog on until it is gone.
END SUBROUTINE BIGADDN !Negative values for B are not considered.
SUBROUTINE BIGMULTN(B,N) !B:=B*N; Multiply by an integer possibly bigger than the base.
TYPE(BIGNUM) B !The worker.
INTEGER N !A computer number, not a multi-digit number.
INTEGER D !Must be able to hold (BIGBASE - 1)*N + C
INTEGER C !The carry to the next digit.
INTEGER*8 DD !N may be large...
INTEGER V !Beware of BIGMULT(B,B); B.DIGIT(1) is N.
INTEGER I !A stepper.
IF (N.EQ.1) RETURN !Phooey.
IF (N.EQ.0) THEN !This takes a little more effort.
B.LAST = 1 !One digit.
B.DIGIT(1) = 0 !Zero.
RETURN !Done.
END IF !Otherwise, start work.
V = N !Grab a local copy in case of doubled reference.
C = 0 !No previous digit to carry from.
IF (N .LT. HUGE(D)/BIGBASE) THEN !Can D hold N*BIGBASE?
DO I = 1,B.LAST !Step through the digits, upwards powers.
D = B.DIGIT(I)*V + C !Grab a digit and apply the multiply.
B.DIGIT(I) = MOD(D,BIGBASE) !Place the resulting digit.
C = D/BIGBASE !Agony! TWO divisions per step!!
END DO !On to the next digit up.
ELSE !Larger N means escalating to a proper double-digit product.
DO I = 1,B.LAST !Step through the same digits.
DD = INT8(B.DIGIT(I))*V + C !Grab a digit and apply the multiply.
B.DIGIT(I) = MOD(DD,BIGBASE) !Place the resulting digit.
C = DD/BIGBASE !Agony! TWO divisions per step!!
END DO !On to the next digit up.
END IF !Either way, the carry is a single digit.
DO WHILE(C .GT. 0) !Now spread the last carry to further digits.
B.LAST = B.LAST + 1 !Up one more.
IF (B.LAST .GT. BIGLIMIT) STOP "Overflow by multiply!" !Perhaps not.
B.DIGIT(B.LAST) = MOD(C,BIGBASE) !The digit.
C = C/BIGBASE !The carry may be large, if N is large.
END DO !So slog on until it is gone.
END SUBROUTINE BIGMULTN !Primary school stuff.
SUBROUTINE BIGMULT(A,B) !A:=A*B, and yes, BIGMULT(A,A) will square A.
C Calculates from the high-order end downwards, with carries going against that flow. C This enables the result to be calculated in-place, even with BIGMULT(A,A). C C a5 a4 a3 a2 a1 ...Five digits. c x b4 b3 b2 b1 ...Four digits. c ----------------------------------------- c a5b1 a4b1 a3b1 a2b1 a1b1 c a5b2 a4b2 a3b2 a2b2 a1b2 ... five wide. c a5b3 a4b3 a3b3 a2b3 a1b3 ... four high. c a5b4 a4b4 a3b4 a2b4 a1b4 c --------------------------------------------- c carry 8 7 6 5 4 3 2 1 At most nine digits. c --------------------------------------------- C A column sum of many digit products might overflow S given many B digits and a big base.
TYPE(BIGNUM) A,B !The numbers.
INTEGER IA,FA,LA !For A: Index, First (highest order), Last (lowest order).
INTEGER IB,FB,LB !For B: Index, First (highest order), Last (lowest order).
INTEGER L !Fingers a digit for a result.
INTEGER*8 S !Scratchpad for digit multiply and summation.
IF ((BIGBASE - 1)*B.LAST.GE.HUGE(S)/(BIGBASE - 1)) !Max. digit product, summed,
1 STOP "BIGMULT: too many B digits! Could overflow S!" !Rather than only ever one at a time..
Check for some simple situations in the hope of evading big nullities.
LB = B.LAST !I'll need a copy of the original layout.
IF (LB.EQ.1) THEN !A single digit B may have some special values.
IF (B.DIGIT(1).EQ.0) THEN !Multiplying A by zero?
CALL BIGLOADN(A,0) !Easy.
ELSE IF (B.DIGIT(1).NE.1) THEN !And if we're not multiplying by one,
CALL BIGMULTN(A,B.DIGIT(1)) !This is easier.
END IF !These values do not appear at random.
RETURN !Done.
END IF !Single-digit big numbers? Hummm.
LA = A.LAST !Copy for later also.
IF (LA.EQ.1) THEN !Perhaps a single-digit here instead..
IF (A.DIGIT(1).EQ.0) RETURN !Zero times something? Zero!
IF (A.DIGIT(1).EQ.1) THEN !One times something?
A.LAST = B.LAST !Yes!
A.DIGIT(1:A.LAST) = B.DIGIT(1:B.LAST) !Copy that something.
RETURN !That was easy.
END IF !So much for one.
END IF !Might as well avoid wasted effort.
Can't avoid work any longer. If the big numbers are *big*, then a lot can be avoided. But, single-digits surely would be rare?
IF (LA + LB .GT. BIGLIMIT + 1) STOP "BigMult will overflow!" !Fixed storage sizes.
A.LAST = LA + LB - 1 !Where the first result digit will go.
FA = LA !Parallelogram parsing control.
FB = LB !These finger the high-order digit.
Commence producing digits. Work down each column, starting with the high-order side.
S = 0 !The grand sum, of many digit products.
10 IB = FB !Index of B's digit of the moment.
DO IA = FA,LA,-1 !Accumulate, working down a column, though bottom to top would do also.
S = INT8(A.DIGIT(IA))*B.DIGIT(IB) + S !Another digit product added in.
IB = IB + 1 !NB: IA + IB is constant in this loop.
END DO !It is the digit power, plus two since DIGIT arrays start with one.
IF (S.LT.0) STOP "BIGMULT: S has flipped sign!" !Oh dear.
L = FA + FB - 1 !Finger the recipient digit.
A.DIGIT(L) = MOD(S,BIGBASE) !Place the digit.
Completed the column sum. Now spread any carry into higher digits as necessary.
S = S/BIGBASE !The sum's carry to the next digit up.
DO WHILE(S > 0) !It may well be multi-digit itself, being the sum of many digit products.
L = L + 1 !Go back (up) a power.
IF (L.LE.A.LAST) THEN !Since we're going high to low,
S = S + A.DIGIT(L) !We add to prior work as we go back up.
ELSE !Or else (once, sigh), extending.
A.LAST = L !Because the high-order product carried up one.
IF (L.GE.BIGLIMIT) STOP "BigMult has overflowed!" !Perhaps too far!
END IF !Righto, S is ready.
A.DIGIT(L) = MOD(S,BIGBASE) !Place the digit.
S = S/BIGBASE !Drop a power as we climb to a still higher digit.
END DO !This may be many powers large.
Contemplate the parallelogram.
IF (FB.GT.1) THEN !The topmost term of a column
FB = FB - 1 !Starts with this B-digit.
ELSE IF (FA.GT.1) THEN !But after the units B-digit is reached,
FA = FA - 1 !The column starts with lesser A-digits.
ELSE !And when the units A-digit has been reached,
RETURN !We're done.
END IF !So much for the start elements of the parallelogram.
IF (LA.GT.1) LA = LA - 1 !The lowest-order A-digit of a column.
GO TO 10 !Peruse the diagram.
END SUBROUTINE BIGMULT !Not the Primary School order, but equivalent.
SUBROUTINE BIGSQUARE(B) !B:=B*B.
c The special feature here is that half the effort of digit multiplying can be avoided by noting c that many digit products are paired, since Bi*Bj = Bj*Bi. Working along the rows from right to left c in the usual manner is easy: start with the Bi*Bi term, then roll along adding 2*Bi*Bj terms in long arithmetic. c The loop control is simple enough, but the doubling requirement is annoying since with fully-packed words c there would be overflows to deal with. Another way is to go down the columns. c "But you still have to double the terms" said Bruce Christianson, whom I met on the way back from lunch. c "Yes, but this way you can add them all up and double them once at the end." c Howl of anguish from Bruce, bewildered blinking from me: "What did I say?" c Bruce had spent most of a lunchtime discussing with John Rumsey, VAX expert, the opportunities offered c by Digital's VAX cpu, especially the feature that allows custom microcode to be loaded for unused op-codes c and which might help with the irritating task of the doubling of each term... c Actually, this can be done in yet a third way, in two passes: first, form the result (with no doubling) c and no squared terms (that are not to be doubled) in any way convenient, then in the second pass, c perform the doubling and add in the squared terms. Running along the rows has the advantage that c one of the digits is constant, so only the other has to be extracted from the array. c But, results must be placed in array elements, so there would still be two array accesses per step. c C C b6 b5 b4 b3 b2 b1 ... Six digits. c x b6 b5 b4 b3 b2 b1 c ----------------------------------------------------------------------- c b6b1 b5b1 b4b1 b3b1 b2b1 b1b1 c b6b2 b5b2 b4b2 b3b2 b2b2 b1b2 ... six wide. c b6b3 b5b3 b4b3 b3b3 b2b3 b1b3 ... six high. c b6b4 b5b4 b4b4 b3b4 b2b4 b1b4 c b6b5 b5b5 b4b5 b3b5 b2b5 b1b5 c b6b6 b5b6 b4b6 b3b6 b2b6 b1b6 c c But, since b2b1 = b1b2, etc., c c |2b6b1 2b5b1|2b4b1 2b3b1|2b2b1 b1b1 ... six, c 2b6b2|2b5b2 2b4b2|2b3b2 b2b2| ... five, c |2b6b3 2b5b3|2b4b3 b3b3| | ... four, c 2b6b4|2b5b4 b4b4| | | ... three, c |2b6b5 b5b5| | | | ... two, c b6b6| | | | | ... one. c ----------------------------------------------------------------------- c carry 11 10 9 8 7 6 5 4 3 2 1 At most twelve digits. c ----------------------------------------------------------------------- C A column sum of many digit products might overflow S given many B digits and a big base. c Working from the high-order end downwards (with carries going against that flow) enables the result to be calculated in-place.
TYPE(BIGNUM) B !The number to be squared.
INTEGER FA,LA !Index for the first and last of the lead digit in the mixed digit products. a6 in a6a5, etc.
INTEGER FB !Index for the first of the second digit in the product pairs.
INTEGER IA,IB !Index for the first and second parts of each digit product a5a4, etc.
INTEGER L !Fingers a digit for a result.
INTEGER*8 S !Scratchpad for digit multiply and summation.
LOGICAL FLIP !Some columns have a lone term, others have only doubled terms.
c write (6,*) "B-digit limit",HUGE(S)/INT8(BIGBASE - 1)**2
IF ((BIGBASE - 1)*B.LAST.GE.HUGE(S)/(BIGBASE - 1)) !Max. digit product, summed,
1 STOP "BIGSQUARE: too many B digits! Could overflow S!" !Rather than only ever one at a time..
Check for some simple situations in the hope of evading big nullities.
FB = B.LAST !I'll need a copy of the original layout.
IF (FB.EQ.1) THEN !A single digit B doesn't have many digits.
CALL BIGMULTN(B,B.DIGIT(1)) !This is easier.
RETURN !Done.
END IF !Single-digit big numbers? Hummm.
Can't avoid work any longer. If the big numbers are *big*, then a lot can be avoided. But, single-digits surely would be rare?
L = 2*FB - 1 !Locate the highest order product: b6b6 in the schedule.
IF (L .GT. BIGLIMIT) STOP "BigSquare will overflow!" !Fixed storage sizes.
B.DIGIT(FB + 1:L) = 0 !Scrub, ready for carry propagation.
B.LAST = L !This might be extended by one.
FA = FB !Triangle parsing control. These finger the high-order digit.
LA = FA + 1 !Syncopation: the a6a6 lone term has no sum of doubled products.
FLIP = .TRUE. !So the first loop for them won't. L is odd.
Commence producing digits. Work down each column, starting with the high-order side.
S = 0 !The grand sum, of many digit products.
10 IB = FB !Index of the second digit in each product.
DO IA = FA,LA,-1 !Accumulate, working down a column, though bottom to top would do also.
S = INT8(B.DIGIT(IA))*B.DIGIT(IB) + S !Another digit product added in.
IB = IB + 1 !NB: IA + IB is constant in this loop.
END DO !It is the digit power, plus two since DIGIT arrays start with one..
S = S + S !Thus the sum of the 2* terms in one go...
IF (FLIP) THEN !For odd powers, there is also an undoubled product.
LA = LA - 1 !Its index is one less than the last digit for a doubled product.
S = S + INT8(B.DIGIT(LA))**2 !Another chance for overflowing S.
END IF !The column sum is now complete.
IF (S.LT.0) STOP "BigSquare: S has flipped sign!" !Oh dear.
FLIP = .NOT.FLIP !Thus the adjustment on every odd power of the result.
L = FA + FB - 1 !Finger the recipient digit.
B.DIGIT(L) = MOD(S,BIGBASE) !Place the digit.
Completed the column sum. Now spread any carry into higher digits as necessary.
S = S/BIGBASE !The sum's carry to the next digit up.
DO WHILE(S > 0) !It may well be multi-digit itself, being the sum of many digit products.
L = L + 1 !Go back (up) a power.
IF (L.LE.B.LAST) THEN !Since we're going high to low,
S = S + B.DIGIT(L) !We add to prior work as we go back up.
ELSE !Or else (once, sigh), extending.
B.LAST = L !Because the high-order product carried up one.
IF (L.GE.BIGLIMIT) STOP "BIGSQUARE has overflowed!" !Perhaps too far.
END IF !Righto, S is ready.
B.DIGIT(L) = MOD(S,BIGBASE) !Place the digit.
S = S/BIGBASE !Drop a power as we climb to a still higher digit.
END DO !This may be many powers large for big summations.
Contemplate the parallelogram.
IF (FB.GT.1) THEN !The topmost term of a column
FB = FB - 1 !Starts with this second digit of each product pair.
ELSE IF (FA.GT.1) THEN !But after the units is reached,
FA = FA - 1 !The column starts with lesser first digits.
ELSE !And when the units first digit has been reached,
RETURN !We're done.
END IF !So much for the start elements of the parallelogram.
GO TO 10 !Peruse the diagram.
END SUBROUTINE BIGSQUARE !Not the Primary School order, but equivalent.
INTEGER FUNCTION BIGMOD2(B) !B mod 2.
TYPE(BIGNUM) B !The number.
Compilers ought to notice that BIGBASE is a constant, so MOD(BIGBASE,2).EQ.0 is also constant, and not generate pointless code... Could also hope that MOD 2 will be effected via an "and" on binary computers.
IF (MOD(BIGBASE,2).EQ.0) THEN !If the base is even,
BIGMOD2 = MOD(B.DIGIT(1),2) !Only the units digit need be considered.
ELSE !But for odd bases, each digit must be inspected.
BIGMOD2 = MOD(COUNT(MOD(B.DIGIT(1:B.LAST),2) .EQ. 1),2) !Whee!
END IF !That was fun!
END FUNCTION BIGMOD2 !This should be swifter than BIGMODN(B,2)
INTEGER FUNCTION BIGMODN(B,N) !Calculate MOD(B,N), for ordinary N.
Causes 32-bit overflow if any C*BIGBASE + DIGIT(i) exceeds the 32-bit integer limit.
TYPE(BIGNUM) B !The big number, presumed positive.
INTEGER N !The divisor, presumed positive.
INTEGER I !The stepper.
INTEGER C !The carry. Values are in 0:N - 1.
IF (N.LE.0) STOP "BIGMODN: positive divisors only!" !Grr.
C IF (N.GT.46340) STOP "BIGMODN: N limited to 46340!" !If INT8, etc. is unavailable.
C = 0 !Here we go. carry + digit
IF (N.GE.HUGE(C)/BIGBASE) THEN !Maximum term is (N - 1)*BIGBASE + (BIGBASE - 1).
DO I = B.LAST,1,-1 !So, N*BIGBASE - 1 must not exceed HUGE(C).
C = MOD(C*INT8(BIGBASE) + B.DIGIT(I),N) !Or, N*BIGBASE .GT. HUGE(C) + 1.
END DO !Next digit down. !Or, N .GT. HUGE(C)/BIGBASE + 1/BIGBASE
ELSE !Thus, if N is not so big, no escalation is needed.
DO I = B.LAST,1,-1 !From the highest-order digit, downwards.
C = MOD(C*BIGBASE + B.DIGIT(I),N) !Whee!
END DO !Next digit down.
END IF !If INTEGER*8 is not available, then a proper BIGMODN is needed.
BIGMODN = C !The remainder.
END FUNCTION BIGMODN !A simple idea, made complex.
SUBROUTINE BIGMOD(B,M) !B:=B mod M.
c Clear? Hah! Why a four year old child could understand this. c Run out and find me a four year old child, I can't make head or tail out of it. C Groucho Marx as Rufus T. Firefly in Duck Soup. 18min.
TYPE(BIGNUM) B !The number, and the result.
TYPE(BIGNUM) M !The modulus.
INTEGER L !Location in B of M's high-order digit.
INTEGER IB,IM !Fingers to the digits of B and M.
INTEGER Q,P !Quotient and provisional reversal.
INTEGER C !Carry.
INTEGER*8 DD,MM!Working numbers. BIGBASE might be large...
IF (HUGE(DD).LT.FLOAT(BIGBASE - 1)**3) STOP "BIGBASE too big!" !Triple digit size.
Could be easy...
IF (M.LAST.LE.1) THEN !Perhaps a single-digit divisor?
B.DIGIT(1) = BIGMODN(B,M.DIGIT(1)) !Yes! This is much easier.
B.LAST = 1 !A single-digit result.
RETURN !Done. Digits above B.LAST should not be referenced.
END IF !Otherwise, the fun begins. By here, know that M.LAST > 1.
MM = M.DIGIT(M.LAST)*INT8(BIGBASE) + M.DIGIT(M.LAST - 1) !Min. is BIGBASE, = 10 in every base.
Could still be easy...
1 IF (BIGSIGN(B,M)) 3,2,10 !Sign of B - M.
2 B.LAST = 1 !B = M.
B.DIGIT(1) = 0 !Unlikely, but the result is clear.
3 RETURN !B < M: no change.
Can't evade the calculation via detection of silly parameters. Now know that B > M, and so B also has at least two digits.
10 L = B.LAST !Align M to the high-order digit of B.
DD = B.DIGIT(L)*INT8(BIGBASE) + B.DIGIT(L - 1) !Combine the top two digits of B. Max. is (b - 1)*b + (b - 1) = b² - 1.
Q = DD/MM !Estimate the quotient. Max. is BIGBASE - 1, because MM is two digits. Max is (max of DD)/(min of MM; b) = b - 1/b.
IF (Q.LE.0) THEN !E.g. in b = 10, 100 mod 11: 10/11 gives Q = 0. Thus know DD < MM.
L = L - 1 !Move M one over. Must be possible, because B > M.
Q = (DD*BIGBASE + B.DIGIT(L - 1))/MM !But this requires three-digit arithmetic.
END IF !A quotient has been chosen.
Calculate B = B - Q*M, with M aligned to the high order of B via L.
20 C = 0 !Clear the carry, which for subtraction is a "borrow".
IB = L - M.LAST + 1 !Finger the lowest-order in B that is aligned with M's lowest order.
DO IM = 1,M.LAST !Step up through the digits of M.
DD = C + B.DIGIT(IB) - INT8(Q)*M.DIGIT(IM) !Subtract Q*M. Q may be BIGBASE - 1 so a double size result.
IF (DD.LT.0) THEN !Is a borrow needed?
C = (DD + 1)/BIGBASE - 1 !Yes. Perhaps many, if Q > 1. C is negative.
B.DIGIT(IB) = DD - C*BIGBASE !Borrowing -C lots of BIGBASE produces one digit.
ELSE !Otherwise, a positive digit.
B.DIGIT(IB) = DD !No need for MOD(DD,BIGBASE).
C = 0 !Clear any carry.
END IF !So much for that digit.
IB = IB + 1 !Advance one digit up in B.
END DO !And also one digit up in M.
IF (IB.EQ.B.LAST) THEN!If M had been shifted one over,
DD = C + B.DIGIT(B.LAST) !The top digit in B has no corresponding Q*M digit
IF (DD.LT.0) THEN !A borrow for the topmost digit?
C = (DD + 1)/BIGBASE - 1 !Yes. This many.
B.DIGIT(IB) = DD - C*BIGBASE !Correct the digit.
ELSE !But without a borrow,
B.DIGIT(IB) = DD !No need for MOD(DD,BIGBASE)
C = 0 !Clear any previous carry.
END IF !So much for the lonely topmost digit.
END IF !The carry may not be zero, and if so, upper digits of B are complemented in BIGBASE.
Check for overflow. Q might have been too big, when using M instead of just MM. Hopefully, not by much. Consider 1005/101 in base 10. DD = 10, MM = 10 so Q = 1 above: proceed. But, subtracting 101 from 100|5 overflows. Could add back the 101 then shift (and try afresh), or, shift and add back. Quotient is down 10 up 1, thus, down 9.
30 IF (C.LT.0) THEN !Has there been a borrow for the highest digit?
31 C = C*INT8(BIGBASE) + B.DIGIT(L) !Yes. Reconcile the two.
IF (C.EQ.-1 .AND. L.GT.M.LAST) THEN !The borrow came from below?
L = L - 1 !Yes. Shift M along one, if not already all the way.
GO TO 31 !And reconcile further. B - Q*M may clear more than one digit.
END IF !C is negative, and digit-sized.
DD = C*INT8(BIGBASE) + B.DIGIT(L - 1) !The overshoot, a negative number.
P = 1 - DD/MM !Convert to multiples of MM, possibly shifted further.
IB = L - M.LAST + 1 !M's units digit is aligned with this B digit.
C = 0 !Restart, unworried by the existing overflow..
DO IM = 1,M.LAST !Step through the digits of M.
DD = C + B.DIGIT(IB) + INT8(P)*M.DIGIT(IM) !Add in P*M.
B.DIGIT(IB) = MOD(DD,BIGBASE) !Place the digit.
C = DD/BIGBASE !Another divide for the carry.
IB = IB + 1 !Step up one digit in B
END DO !And also one digit in M.
DO IB = IB,B.LAST !Convert any higher inverted digits.
C = C + B.DIGIT(IB) !C is one, and an inverted digit = (BIGBASE - 1).
B.DIGIT(IB) = MOD(C,BIGBASE) !So, this will be zero.
C = C/BIGBASE !And C will be one again.
END DO !All the way to the top of the number.
END IF !So much for any overflow.
Contemplate the result.
CALL BIGNORM(B) !Hopefully, many leading zero digits have appeared...
IF (B.LAST.GT.M.LAST) GO TO 10!If B has more digits, there is definitely more to do.
GO TO 1 !Possibly, B <= M.
END SUBROUTINE BIGMOD !Memories of a hand-cranked adding machine.
INTEGER FUNCTION BIGDIVRN(B,N) !B:=B/N; but, returns the remainder too.
TYPE(BIGNUM) B !The big number, presumed positive, to be divided by N..
INTEGER N !The divisor, presumed positive.
INTEGER I !The stepper.
INTEGER R !The remainder. Values are in 0:N - 1.
INTEGER D !A double digit, if N is not too big..
INTEGER*8 DD !For use with the larger values of N.
c IF (N.GT.46340) STOP "BIGDIVRN: N limited to 46340!" !If INT8, etc. is unavailable.
R = 0 !Here we go.
IF (N.GE.HUGE(D)/BIGBASE) THEN!Edging towards multi-precision arithmetic?
DO I = B.LAST,1,-1 !Step down the digits.
DD = R*INT8(BIGBASE) + B.DIGIT(I) !Form the current digit.
B.DIGIT(I) = DD/N !Place it.
R = MOD(DD,N) !Remainder = R - B.DIGIT(I)*N, in primary school.
END DO !On to the next digit.
ELSE !If a remainder*BIGBASE + digit can fit into D, no need for DD.
DO I = B.LAST,1,-1 !Step down the digits.
D = R*BIGBASE + B.DIGIT(I)!Form the current digit.
B.DIGIT(I) = D/N !Place it.
R = MOD(D,N) !Remainder = D - B.DIGIT(I)*N, in primary school.
END DO !On to the next digit.
END IF !In larger or smaller gulps.
CALL BIGNORM(B) !Trim any leading zero digits.
BIGDIVRN = R !Pass back the remainder as well.
END FUNCTION BIGDIVRN !A pity that two divides per step can't be dodged outside assembler.
SUBROUTINE BIGDIVN(B,N) !B:=B div N !Divide a big number by an ordinary number.
TYPE(BIGNUM) B !The big number to divide.
INTEGER N !The divisor, an ordinary number.
INTEGER D,R !Scratchpads for the calculation.
INTEGER*8 DD !A double-size digit is available...
INTEGER I !A stepper.
IF (N.LE.0) STOP "BIGDIVN: positive divisors only!" !Grr.
C IF (N.GT.46340) STOP "BIGDIVN: N limited to 46340!" !If INT8, etc. is unavailable.
IF (N.EQ.1) GO TO 10 !Nothing need be done, but normalisation might be good.
R = 0 !No higher-order remainder.
IF (N.GE.HUGE(D)/BIGBASE) THEN!Edging towards multi-precision arithmetic?
DO I = B.LAST,1,-1 !Starting at the high-order end and working down.
DD = R*INT8(BIGBASE) + B.DIGIT(I) !A term, double-digit size.
B.DIGIT(I) = DD/N !Divide, losing any remainder.
R = MOD(DD,N) !Max. remainder is N - 1.
END DO !On to the next digit, multiplying the remainder by BIGBASE.
ELSE !With N not large, (N - 1)*BIGBASE + DIGIT won't overflow.
DO I = B.LAST,1,-1 !Starting at the high-order end and working down.
D = R*BIGBASE + B.DIGIT(I) !A term, double-digit size being not too big for D.
B.DIGIT(I) = D/N !Divide, losing any remainder.
R = MOD(D,N) !Recover the remainder.
END DO !On to the next digit.
END IF !Either way, this is as taught in primary school, except there R = D - B.DIGIT(I)*N.
10 CALL BIGNORM(B) !Check for leading zero digits, since B will have been reduced.
END SUBROUTINE BIGDIVN !The remainder could be returned as well.
INTEGER FUNCTION MODEXP(N,X,P) !Calculate X**P mod N without overflowing...
C Relies on a.b mod n = (a mod n)(b mod n) mod n
INTEGER N,X,P !All presumed positive, and X < N.
INTEGER I !A stepper.
INTEGER*8 V,W !Broad scratchpads, otherwise N > 46340 may incur overflow in 32-bit.
V = 1 !=X**0
IF (P.GT.0) THEN !Something to do?
I = P !Yes. Get a copy I can mess with.
W = X !=X**1, X**2, X**4, X**8, ... except, all are mod N.
1 IF (MOD(I,2).EQ.1) V = MOD(V*W,N) !Incorporate W if the low-end calls for it.
I = I/2 !Used. Shift the next one down.
IF (I.GT.0) THEN !Still something to do?
W = MOD(W**2,N) !Yes. Square W ready for the next bit up.
GO TO 1 !Consider it.
END IF !Don't square W if nothing remains. It might overflow.
END IF !Negative powers are ignored.
MODEXP = V !Done, in lb(P) iterations!
END FUNCTION MODEXP !"Bit" presence by arithmetic: works for non-binary arithmetic too.
SUBROUTINE BIGMODEXP(V,N,X,P) !Calculate V = X**P mod N without overflowing...
C Relies on a.b mod n = (a mod n)(b mod n) mod n
TYPE(BIGNUM) V,N,X,P !All presumed positive.
TYPE(BIGNUM) I !A stepper.
TYPE(BIGNUM) W !Broad scratchpads, otherwise N > 46340 may incur overflow in 32-bit.
CALL BIGLOADN(V,1) !=X**0
IF (P.LAST.GT.1 .OR. P.DIGIT(1).GT.0) THEN !Something to do?
I.LAST = P.LAST !Yes. Get a copy I can mess with.
I.DIGIT(1:I.LAST) = P.DIGIT(1:P.LAST) !Only copying the digits in use.
W.LAST = X.LAST !=X**1, X**2, X**4, X**8, ... except, all are mod N.
W.DIGIT(1:W.LAST) = X.DIGIT(1:X.LAST) !Used according to the bits in P.
1 IF (BIGMOD2(I).EQ.1) THEN !Incorporate W if the low-end calls for it.
CALL BIGMULT(V,W) !V:=V*W ...
CALL BIGMOD(V,N) ! ... mod N.
END IF !So much for that bit.
CALL BIGDIVN(I,2) !Used. Shift the next one down.
IF (I.LAST.GT.1 .OR. I.DIGIT(1).GT.0) THEN !Still something to do?
CALL BIGSQUARE(W) !Yes. Square W ready for the next bit up.
CALL BIGMOD(W,N) !Reduced modulo N.
GO TO 1 !Consider it.
END IF !Don't square W if nothing remains. A waste of effort.
END IF !Negative powers are ignored.
END SUBROUTINE BIGMODEXP !"Bit" presence by arithmetic: works for non-binary arithmetic too.
LOGICAL FUNCTION BIGMRPRIME(N,TRIALS) !Could N be a prime number?
USE DFPORT !To get RAND, which returns a "random" number in 0.0 to 1.0. Inclusive.
TYPE(BIGNUM) N !The number to taste.
INTEGER TRIALS !The count of trials to make.
INTEGER S !Counts powers of two in N - 1.
TYPE(BIGNUM) NL1 !Holds N - 1 for comparisons.
TYPE(BIGNUM) D !N - 1 with twos divided out.
TYPE(BIGNUM) A !A number in [2:N - 1]. Any number...
TYPE(BIGNUM) X !Scratchpad.
INTEGER A1,N1 !Assistants for juggling a "random" number.
INTEGER TRIAL,R !Counters.
IF (BIGBASE.LE.3) STOP "BIGMRPRIME: BigBase too small!" !Multi-digit even for small numbers!
Catch some annoying cases.
IF (N.LAST.EQ.1) THEN !A single-digit number?
IF (N.DIGIT(1).LE.4) THEN !Yes. Some special values are known.
BIGMRPRIME = N.DIGIT(1).GE.2 .AND. N.DIGIT(1).LE.3 !Like, the neighbours.
RETURN !Thus allow 2 to be reported as prime.
END IF !Yet, test for 2 as a possible factor for larger numbers.
END IF !Without struggling over SQRT and suchlike.
BIGMRPRIME = .FALSE. !Most numbers are not primes.
IF (BIGMOD2(N).EQ.0) RETURN !A single expression using .OR. risks always evaluating BOTH parts, damnit,
IF (BIGMODN(N,3).EQ.0) RETURN !Even for even numbers. Possibly doing so "in parallel" is no consolation.
Construct D such that N - 1 = D*2**S. By here, N is odd, and greater than three.
D.LAST = N.LAST !Could just put D = N,
D.DIGIT(1:D.LAST) = N.DIGIT(1:D.LAST) !But this copies only the digits in use.
CALL BIGADDN(D,-1) !Thus, D becomes an even number.
NL1.LAST = D.LAST !For later testing of X against N - 1,
NL1.DIGIT(1:NL1.LAST) = D.DIGIT(1:D.LAST) !Retain N - 1.
N1 = MIN(20000000D0,BIGVALUE(NL1)) !Maximum value is N - 1.
S = 1 !Since D is even, it has at least one power of two.
10 CALL BIGDIVN(D,2) !Divide out a power of two.
IF (BIGMOD2(D).EQ.0) THEN !If there is another,
S = S + 1 !Count it,
GO TO 10 !And divide it out also.
END IF !So, D is no longer even. N - 1 = D*2**S
Convince through repetition...
T:DO TRIAL = 1,TRIALS !Some trials come to a definite result.
BIGMRCOUNT(TRIAL) = BIGMRCOUNT(TRIAL) + 1 !
A1 = RAND(0)*(N1 - 2) + 2 !For small N, the birthday problem. NB! RAND can generate 1.
CALL BIGLOADN(A,A1) !A1 is in (0 + 2) = 2 to N - 1 = (1*(N1 - 2) + 2).
CALL BIGMODEXP(X,N,A,D) !X = A**D mod N.
IF (X.LAST.EQ.1 .AND. X.DIGIT(1).EQ.1) CYCLE T !Pox. A prime yields these: 1 or NL1.
IF (BIGSIGN(X,NL1).EQ.0) CYCLE T !A test with .OR. might always evaluate both, damnit.
DO R = 1,S - 1 !Step through the powers of two in N - 1.
CALL BIGSQUARE(X) !X**2 ...
CALL BIGMOD(X,N) ! ... mod N
IF (X.LAST.EQ.1 .AND. X.DIGIT(1).EQ.1) RETURN !X = 1? Definitely composite. No prime does this.
IF (BIGSIGN(X,NL1).EQ.0) CYCLE T !Pox. Try something else.
END DO !Another power of two?
RETURN !Definitely composite.
END DO T !Have another try.
BIGMRPRIME = .TRUE. !Would further trials yield greater assurance?
END FUNCTION BIGMRPRIME !Are some numbers resistant to this scheme?
INTEGER FUNCTION BIGFACTOR(B,LIMIT) !Find the first factor of B. Limited assault.
Careful. The 32-bit integer limit is 2,147,483,647. P(1358124) = 21,474,829 and BIGBASE might be 100.
USE PRIMEBAG !I want a supply.
TYPE(BIGNUM) B !The number to factorise.
INTEGER LIMIT !A limit on the search.
REAL*8 S !The square root of B.
INTEGER F !A candidate factor.
S = SQRT(BIGVALUE(B)) !Pious hope for accuracy...
F = 2 !Here we go.
DO WHILE(F.LE.S) !Thus avoid declaring 2 to have a factor of 2.
IF (BIGMODN(B,F).EQ.0) THEN !Divisible?
BIGFACTOR = F !Yes!
RETURN !Done.
END IF !But otherwise,
F = NEXTPRIME(F) !Find the next prime number.
IF (F.LE.0 .OR. F.GT.LIMIT) THEN !Overflow? Or, too far already?
BIGFACTOR = 0 !Alas. Sez: Don't know.
c WRITE (6,1) LIMIT,F,S !Make a report.
1 FORMAT ("Passed the limit of ",I0," with ",I0,
1 ", but Sqrt(n) =",E16.8)
IF (BIGMRPRIME(B,BIGMRTRIALS)) THEN !Crank up a Miller-Rabin probe.
BIGFACTOR = -1 !Non-priminess not observed in this probe.
ELSE !Or possibly,
BIGFACTOR = -2 !Non-primeness definitely observed.
END IF !That was mysterious.
RETURN !Pass the word..
END IF !But otherwise,
END DO !Consider the next factor.
BIGFACTOR = 1 !If we pass SQRT(B), B is definitely prime.
END FUNCTION BIGFACTOR !The first factor suffices.
END MODULE BIGNUMBERS !No fancy tricks.
</lang>
Use
With all that in hand, big numbers can be manipulated, and function BIGFACTOR can be used where others use something like IsProbablyPrime. So the task is solvable via <lang Fortran> PROGRAM PRIMORIALP !Simple enough, with some assistants.
USE PRIMEBAG !Some prime numbers are wanted.
USE BIGNUMBERS !Just so.
TYPE(BIGNUM) B !P'll have one.
INTEGER MAXF !Largest factor to consider by direct division.
PARAMETER (MAXF = 1000000) !Some determination.
INTEGER I,L,P !Step stuff.
INTEGER FU,FD,FUD(2) !Found factors.
EQUIVALENCE (FUD(1),FD),(FUD(2),FU) !A struggle to prepare an array.
INTEGER NHIT,HIT(666) !A little list.
CHARACTER*4 WOT !A remark.
CHARACTER*66 ALINE !A scratchpad.
REAL T0,T1 !In memory of lost time.
MSG = 6 !Standard output.
WRITE (MSG,1) BIGLIMIT,BIGBASE,HUGE(P),HUGE(INT8(P)),MAXF !Announce.
1 FORMAT ('Calculates primorial "primes"',/,
1 "A primorial prime is a value N such that",/,
2 " Primorial(N) - 1 is prime, OR",/,
3 " Primorial(N) + 1 is prime, or both.",/,
4 "and Primorial(N) is the product of the first N prime numbers.",/
5 "Working with up to ",I0," digits in base ",I0,"."/
6 "Integer limits are ",I0," and ",I0/
7 "Primes up to ",I0," are tried as possible factors.")
IF (.NOT.GRASPPRIMEBAG(66)) STOP "Gan't grab my file!" !Attempt in hope.
WRITE (MSG,2)
2 FORMAT (34X,"First First"/,
1 "Primorial#",5X,"Approx.",7X," -1 Factor +1 Factor Hit")
Commence prime mashing.
100 NHIT = 0 !My list is empty.
B.LAST = 1 !Begin at the beginning.
B.DIGIT(1) = 1 !With one. The same, whatever BIGBASE.
CALL CPU_TIME(T0) !Start the timing.
DO P = 1,80 !460 !Step along the primorials. 457 is the 20'th.
CALL BIGMULTN(B,PRIME(P)) !Multiply by the next prime.
c WRITE (MSG,101) P,PRIME(P),P,B.DIGIT(B.LAST:1:-1) !Digits in Arabic/Hindu order.
101 FORMAT ("Prime(",I0,") = ",I0,", Primorial(",I0,") = ", !For a possibly multi-BIGBASE sequence.
1 I0,29I<BIGORDER>.<BIGORDER>,/,(10I<BIGORDER>.<BIGORDER>)) !The first without leading zero digits.
FU = 0 !No factor for up one.
FD = 0 !No factor for down one.
CALL BIGADDN(B,+1) !Go up one.
FU = BIGFACTOR(B,MAXF) !Find a factor, maybe.
CALL BIGADDN(B,-2) !Now test down one.
IF (ABS(FU).NE.1) FD = BIGFACTOR(B,MAXF) !But only if FU didn't report "prime".
IF (ABS(FU).EQ.1 .OR. ABS(FD).EQ.1) THEN !Since if either candidate is a prime,
WOT = "Yes!" !Then a hit.
NHIT = NHIT + 1 !So count up a success.
HIT(NHIT) = P !And append to my list.
ELSE IF (ABS(FU).GT.1 .AND. ABS(FD).GT.1) THEN !But if both have factors,
WOT = "No." !Then definitely not a hit.
ELSE !Otherwise,
WOT = "?" !I can't decide.
END IF !So much for that candidate.
CALL BIGADDN(B,1) !Recover the original primorial value.
CALL BIGTASTE(B) !Prepare an approximate value in BIGLEAD.
WRITE (ALINE,102) P,BIGLEAD,FD,FU,WOT !Prepare a report.
102 FORMAT (I10," 0.",<BIGLEADN>I0,T24,"E+",I0,T30,I10,I10,1X,A) !A table.
DO I = 1,2 !Two columns to the table.
L = 20 + I*10 !Each with the same interpretation.
SELECT CASE(FUD(I)) !So, one set of code.
CASE(-2); ALINE(L:L+9) = " Composite"
CASE(-1); ALINE(L:L+9) = " MR prime"
CASE( 0); ALINE(L+9:L+9) = "?" !Don't know. Didn't look.
CASE DEFAULT !All other values stay. The first factor.
END SELECT !So much for re-interpretation.
END DO !On to the next column.
WRITE (MSG,"(A)") ALINE !Show the report.
END DO !On to the next prime.
CALL CPU_TIME(T1) !Completed the run.
Cast forth some pearls.
200 WRITE (MSG,201) !First, some statistics.
201 FORMAT (/,"The MR prime test makes a series of trials, "
1 "stopping early",/'only when a "definitely composite" ',
2 "result is encountered.")
WRITE (MSG,202) "Trial",(I,I = 1,BIGMRTRIALS) !Roll the trial number.
WRITE (MSG,202) "Count",BIGMRCOUNT !Now the counts.
202 FORMAT (A6,": ",666I6) !This should do.
WRITE (MSG,203) NHIT,HIT(1:NHIT) !The list of primorial "primes"..
203 FORMAT (/,I0," in the hit list: ",I0,666(",",I0:)) !Don't actually expect so many.
WRITE (MSG,*) "CPU time:",T1 - T0 !The cost.
END !So much for that.
</lang> A slightly more complex scheme to interpret the more complex mumbling from BIGFACTOR, which is less sure of its reports in more ways.
Output
A somewhat lesser upper bound for factoring leads to a faster run.
Calculates primorial "primes"
A primorial prime is a value N such that
Primorial(N) - 1 is prime, OR
Primorial(N) + 1 is prime, or both.
and Primorial(N) is the product of the first N prime numbers.
Working with up to 2222 digits in base 10000.
Integer limits are 2147483647 and 9223372036854775807
Primes up to 1000000 are tried as possible factors.
First First
Primorial# Approx. -1 Factor +1 Factor Hit
1 0.20 E+1 ? 1 Yes!
2 0.60 E+1 ? 1 Yes!
3 0.300 E+2 ? 1 Yes!
4 0.2100 E+3 ? 1 Yes!
5 0.23100 E+4 ? 1 Yes!
6 0.330 E+5 1 59 Yes!
7 0.51510 E+6 61 19 No.
8 0.9699690 E+7 53 347 No.
9 0.22309 E+9 37 317 No.
10 0.646969 E+10 79 331 No.
11 0.20056049 E+12 ? 1 Yes!
12 0.74207 E+13 229 181 No.
13 0.3042503 E+15 MR prime 61 Yes!
14 0.13083 E+17 141269 167 No.
15 0.614890 E+18 191 953 No.
16 0.32589158 E+20 87337 73 No.
17 0.192276 E+22 Composite 277 No.
18 0.11728838 E+24 1193 223 No.
19 0.78583 E+25 163 Composite No.
20 0.5579408 E+27 Composite 1063 No.
21 0.4730 E+29 313 2521 No.
22 0.3217645 E+31 163 22093 No.
23 0.26706 E+33 139 265739 No.
24 0.2376874 E+35 MR prime 131 Yes!
25 0.23056 E+37 66683 Composite No.
26 0.2328624 E+39 Composite 960703 No.
27 0.23985 E+41 15649 2297 No.
28 0.2566376 E+43 Composite 149 No.
29 0.27973 E+45 719 334507 No.
30 0.3161005 E+47 295201 Composite No.
31 0.4145 E+49 Composite 1543 No.
32 0.5258965 E+51 Composite 1951 No.
33 0.72048 E+53 Composite 881 No.
34 0.10014647 E+56 4871 Composite No.
35 0.149218 E+58 673 Composite No.
36 0.22531953 E+60 311 Composite No.
37 0.353752 E+62 1409 1381 No.
38 0.57661522 E+64 1291 1361 No.
39 0.962947 E+66 331 Composite No.
40 0.16659 E+69 Composite Composite No.
41 0.2981959 E+71 23497 601 No.
42 0.53973 E+73 711427 107453 No.
43 0.10308931 E+76 521 32999 No.
44 0.198962 E+78 673 Composite No.
45 0.39195588 E+80 519577 521831 No.
46 0.779992 E+82 Composite 467 No.
47 0.16458 E+85 56543 1051 No.
48 0.36797 E+87 811 Composite No.
49 0.83311 E+89 182309 3187 No.
50 0.19078267 E+92 53077 126173 No.
51 0.444524 E+94 641 5981 No.
52 0.1624 E+97 349 313 No.
53 0.256412 E+99 389 3499 No.
54 0.64266 E+101 565921 Composite No.
55 0.16516447 E+104 777041 34259 No.
56 0.434383 E+106 Composite 271 No.
57 0.11685 E+109 757 17291 No.
58 0.3166605 E+111 2341 15559 No.
59 0.87715 E+113 44753 243233 No.
60 0.24647906 E+116 Composite 62131 No.
61 0.697536 E+118 7027 307 No.
62 0.2438 E+121 4111 Composite No.
63 0.6274404 E+123 1571 Composite No.
64 0.195134 E+126 Composite Composite No.
65 0.61076929 E+128 Composite 3761 No.
66 0.1936139 E+131 MR prime Composite Yes!
67 0.64086 E+133 Composite 8089 No.
68 0.21597046 E+136 MR prime 38453 Yes!
69 0.749417 E+138 1609 1579 No.
70 0.26155 E+141 33409 44953 No.
71 0.9232599 E+143 Composite 2143 No.
72 0.331450 E+146 1069 Composite No.
73 0.12164 E+149 Composite Composite No.
74 0.4537256 E+151 Composite 4973 No.
75 0.171962 E+154 ? MR prime Yes!
76 0.65861450 E+156 5407 Composite No.
77 0.2562010 E+159 104707 216401 No.
78 0.101712 E+162 599 8629 No.
79 0.40786437 E+164 4139 356647 No.
80 0.1668165 E+167 35447 Composite No.
The MR prime test makes a series of trials, stopping early
only when a "definitely composite" result is encountered.
Trial: 1 2 3 4 5 6
Count: 41 5 5 5 5 5
12 in the hit list: 1,2,3,4,5,6,11,13,24,66,68,75
CPU time: 4.125000
Function BIGFACTOR calls upon BIGMRPRIME to make six trials. Somewhere I saw a note that a MR test might fail to denounce a "random" composite number about a quarter of the time, so, up to six trials perhaps reduces this to less than one in a thousand. Except that none of the tested numbers can be considered as a "random" number! They have special properties! Being one off from a number with a sequence of prime factors, a primorial number, and note that none of the values tested turn out to have a factor of two, or three, or five, when "on average", surely half the "random" numbers would be divisible by two, a third by three, etc. But, being one away from Primorial N, a number which is a multiple of two, three, five, etc. means that they can't be divisible by two, three, five, etc. until a possible factor exceeds N.
Evidently, a candidate is denounced as definitely composite in the first trial, at least for these candidates. The following trials seem wasted effort...
Different-size bases have an effect, though not a regular one.
Base 10 100 1000 10000 100000 1000000
Time 24.55 8.20 4.45 4.11 3.48 2.92 Using 64-bit integers.
11.17 4.09 2.55 Using 32-bit integers.
For software compatibility reasons too irritating to recount, this system runs 32-bit Windows XP even though its AMD Six-Core processor is 64-bit (and when running GNU-Linux, uses 64-bit) thus the INTEGER*8 variables used to handle double-digit products (and even triple-digit products) have their arithmetic done via software, not hardware. While it is amusing to have software for multi-precision integer arithmetic supported via software for double-precision integer arithmetic, this comes at a cost! There seems to be no attempt to use the 64-bit facilities of the floating-point arithmetic unit.
Converting from BIGDIVN(B,2) to BIGDIV2(B) had very little effect on the speed, but inspection of the code demonstrated that the Compaq compiler carried forward the result of expressions involving constants such as IF (MOD(BIGBASE,2).EQ.0) THEN at compile time so that no run-time test was made, and "dead code" for the ELSE option was not generated. Further, MOD(N,2) generated code involving an and, so there is no need to mess about by coding an AND or fiddling with special syntax that would only work on a binary computer. With such a compiler, a sort of "conditional compilation" is available without needing the often grotesque syntax of systems explicitly offering such facilities.
Replacing BIGMULT(B,B) by BIGSQUARE(B) converted ~4·28 seconds to ~4·14 seconds for the test run; three percent.
FreeBASIC
<lang freebasic>' version 23-10-2016 ' compile with: fbc -s console
- Define max 9999 ' max number for the sieve
- Include Once "gmp.bi"
Dim As mpz_ptr p, p1 p = Allocate(Len(__mpz_struct)) : Mpz_init_set_ui(p, 1) p1 = Allocate(Len(__mpz_struct)) : Mpz_init(p1)
Dim As UInteger i, n, x Dim As Byte prime(max)
' Sieve of Eratosthenes For i = 4 To max Step 2
prime(i) = 1
Next For i = 3 To Sqr(max) Step 2
If prime(i) = 1 Then Continue For For n = i * i To max Step i * 2 prime(n) = 1 Next
Next
n = 0 : x = 0 For i = 2 To max
If prime(i) = 1 Then Continue For x = x + 1 mpz_mul_ui(p, p, i) mpz_sub_ui(p1, p, 1) If mpz_probab_prime_p(p1, 25) > 0 Then Print Using "####"; x; : Print ","; n += 1 If n >= 20 Then Exit For Continue For End If mpz_add_ui(p1, p, 1) If mpz_probab_prime_p(p1, 25) > 0 Then Print Using "####"; x; : Print ","; n += 1 If n >= 20 Then Exit For End If
Next
Print mpz_clear(p) mpz_clear(p1)
' empty keyboard buffer While Inkey <> "" : Wend Print : Print "hit any key to end program" Sleep End</lang>
- Output:
1, 2, 3, 4, 5, 6, 11, 13, 24, 66, 68, 75, 167, 171, 172, 287, 310, 352, 384, 457,
Haskell
<lang Haskell> import Data.List (scanl1, findIndices, nub)
primes :: [Integer] primes = 2 : filter isPrime [3,5..]
isPrime :: Integer -> Bool isPrime = isPrime_ primes
where isPrime_ :: [Integer] -> Integer -> Bool
isPrime_ (p:ps) n
| p * p > n = True
| n `mod` p == 0 = False
| otherwise = isPrime_ ps n
primorials :: [Integer] primorials = 1 : scanl1 (*) primes
primorialsPlusMinusOne :: [Integer] primorialsPlusMinusOne = concatMap (\n -> [n-1, n+1]) primorials
sequenceOfPrimorialPrimes :: [Int] sequenceOfPrimorialPrimes = tail $ nub $ map (`div` 2) $ findIndices (==True) bools
where bools = map isPrime primorialsPlusMinusOne
main :: IO () main = mapM_ print $ take 10 sequenceOfPrimorialPrimes </lang>
- Output:
1 2 3 4 5 6 11 13 24 66
J
<lang J>primoprim=: [: I. [: +./ 1 p: (1,_1) +/ */\@:p:@i.</lang>
This isn't particularly fast - J's current extended precision number implementation favors portability over speed.
Example use:
<lang J> primoprim 600x 0 1 2 3 4 5 10 12 23 65 67 74 166 170 171 286 309 351 383 456 563 589</lang>
Note that the task specifies that index zero is not a part of the sequence for this task. So pretend you didn't see it.
Java
<lang java>import java.math.BigInteger;
public class PrimorialPrimes {
final static int sieveLimit = 1550_000; static boolean[] notPrime = sieve(sieveLimit);
public static void main(String[] args) {
int count = 0;
for (int i = 1; i < 1000_000 && count < 20; i++) {
BigInteger b = primorial(i);
if (b.add(BigInteger.ONE).isProbablePrime(1)
|| b.subtract(BigInteger.ONE).isProbablePrime(1)) {
System.out.printf("%d ", i);
count++;
}
}
}
static BigInteger primorial(int n) {
if (n == 0)
return BigInteger.ONE;
BigInteger result = BigInteger.ONE;
for (int i = 0; i < sieveLimit && n > 0; i++) {
if (notPrime[i])
continue;
result = result.multiply(BigInteger.valueOf(i));
n--;
}
return result;
}
public static boolean[] sieve(int limit) {
boolean[] composite = new boolean[limit];
composite[0] = composite[1] = true;
int max = (int) Math.sqrt(limit);
for (int n = 2; n <= max; n++) {
if (!composite[n]) {
for (int k = n * n; k < limit; k += n) {
composite[k] = true;
}
}
}
return composite;
}
}</lang>
1 2 3 4 5 6 11 13 24 66 68 75 167 171 172 287 310 352 384 457
Kotlin
<lang scala>// version 1.0.6
import java.math.BigInteger
const val LIMIT = 20 // expect a run time of about 2 minutes on a typical laptop
fun isPrime(n: Int): Boolean {
if (n < 2) return false
if (n % 2 == 0) return n == 2
if (n % 3 == 0) return n == 3
var d : Int = 5
while (d * d <= n) {
if (n % d == 0) return false
d += 2
if (n % d == 0) return false
d += 4
}
return true
}
fun main(args: Array<String>) {
println("The first $LIMIT primorial indices in the sequence are:")
print("1 ")
var primorial = 1
var count = 1
var p = 3
var prod = BigInteger.valueOf(2L)
while(true) {
if (isPrime(p)) {
prod *= BigInteger.valueOf(p.toLong())
primorial++
if ((prod + BigInteger.ONE).isProbablePrime(1) || (prod - BigInteger.ONE).isProbablePrime(1)) {
print("$primorial ")
count++
if (count == LIMIT) {
println()
break
}
}
}
p += 2
}
}</lang>
- Output:
The first 20 primorial indices in the sequence are: 1 2 3 4 5 6 11 13 24 66 68 75 167 171 172 287 310 352 384 457
PARI/GP
This program, in principle, generates arbitrarily many primirial prime indices. Practically speaking it depends on your patience; it generates them up to 643 in about 5 seconds. <lang parigp>n=0; P=1; forprime(p=2,, P*=p; n++; if(ispseudoprime(P+1) || ispseudoprime(P-1), print1(n", ")))</lang>
- Output:
1, 2, 3, 4, 5, 6, 11, 13, 24, 66, 68, 75, 167, 171, 172, 287, 310, 352, 384, 457, 564, 590, 616, 620, 643, 849,
Perl
<lang perl>use ntheory ":all"; my $i = 0; for (1..1e6) {
my $n = pn_primorial($_);
if (is_prime($n-1) || is_prime($n+1)) {
print "$_\n";
last if ++$i >= 20;
}
}</lang>
- Output:
1 2 3 4 5 6 11 13 24 66 68 75 167 171 172 287 310 352 384 457
Perl 6
<lang perl6>constant @primes = |grep *.is-prime, 2..*; constant @primorials = [\*] 1, @primes;
my @pp_indexes := |@primorials.pairs.map: {
.key if ( .value + any(1, -1) ).is-prime
};
say ~ @pp_indexes[ 0 ^.. 20 ]; # Skipping bogus first element.</lang>
- Output:
1 2 3 4 5 6 11 13 24 66 68 75 167 171 172 287 310 352 384 457
Python
Uses the pure python library pyprimes . Note that pyprimes.isprime switches to a fast and good probabilistic algorithm for larger integers.
<lang python>import pyprimes
def primorial_prime(_pmax=500):
isprime = pyprimes.isprime
n, primo = 0, 1
for prime in pyprimes.nprimes(_pmax):
n, primo = n+1, primo * prime
if isprime(primo-1) or isprime(primo+1):
yield n
if __name__ == '__main__':
# Turn off warning on use of probabilistic formula for prime test
pyprimes.warn_probably = False
for i, n in zip(range(20), primorial_prime()):
print('Primorial prime %2i at primorial index: %3i' % (i+1, n))</lang>
- Output:
Primorial prime 1 at primorial index: 1 Primorial prime 2 at primorial index: 2 Primorial prime 3 at primorial index: 3 Primorial prime 4 at primorial index: 4 Primorial prime 5 at primorial index: 5 Primorial prime 6 at primorial index: 6 Primorial prime 7 at primorial index: 11 Primorial prime 8 at primorial index: 13 Primorial prime 9 at primorial index: 24 Primorial prime 10 at primorial index: 66 Primorial prime 11 at primorial index: 68 Primorial prime 12 at primorial index: 75 Primorial prime 13 at primorial index: 167 Primorial prime 14 at primorial index: 171 Primorial prime 15 at primorial index: 172 Primorial prime 16 at primorial index: 287 Primorial prime 17 at primorial index: 310 Primorial prime 18 at primorial index: 352 Primorial prime 19 at primorial index: 384 Primorial prime 20 at primorial index: 457
Racket
We use a memorized version of primordial to make it faster.
<lang Racket>#lang racket
(require math/number-theory
racket/generator)
(define-syntax-rule (define/cache (name arg) body ...)
(begin
(define cache (make-hash))
(define (name arg)
(hash-ref! cache arg (lambda () body ...)))))
(define/cache (primorial n)
(if (zero? n)
1
(* (nth-prime (sub1 n))
(primorial (sub1 n)))))
(for ([i (in-range 20)]
[n (in-generator (for ([i (in-naturals 1)])
(define pr (primorial i))
(when (or (prime? (add1 pr)) (prime? (sub1 pr)))
(yield i))))])
(displayln n))</lang>
- Output:
1 2 3 4 5 6 11 13 24 66 68 75 167 171 172 287 310 352 384 457
Sidef
<lang ruby>func primorial_primes(n) {
var k = 1 var p = 2 var P = 2
var seq = []
for (var i = 0; i < n; ++k) {
if (is_prime(P-1) || is_prime(P+1)) {
seq << k
++i
}
p.next_prime!
P *= p
}
return seq
}
say primorial_primes(20)</lang>
- Output:
[1, 2, 3, 4, 5, 6, 11, 13, 24, 66, 68, 75, 167, 171, 172, 287, 310, 352, 384, 457]
zkl
Uses libGMP (GNU MP Bignum Library) <lang zkl>var [const] BN=Import("zklBigNum"); // libGMP p,s,i,n:=BN(1),BN(1), 0,0; do{ n+=1;
s.nextPrime(); // in place, probabilistic
p.mul(s); // in place
if((p+1).probablyPrime() or (p-1).probablyPrime()){
println("%3d %5d digits".fmt(n,p.len()));
i+=1;
}
}while(i<20);</lang>
- Output:
1 1 digits 2 1 digits 3 2 digits 4 3 digits 5 4 digits 6 5 digits 11 12 digits 13 15 digits 24 35 digits 66 131 digits 68 136 digits 75 154 digits 167 413 digits 171 425 digits 172 428 digits 287 790 digits 310 866 digits 352 1007 digits 384 1116 digits 457 1368 digits