mathematics

Dec_01_2008 Stumbled across Pascal's triangle again, when multiplying by a scaling matrix. Probably not surprising, but I want to look at it later. >> [1.1*1.1 0 0 0;0 1.1 0 0;0 0 1 0;0 0 0 1] * [1 1 1 1]' ans = 1.2100 1.1000 1.0000 1.0000 >> [11*11 0 0 0;0 11 0 0;0 0 1 0;0 0 0 1] * [1 1 1 1]' ans = 121 11 1 1 Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#PascalsTriangleFromScaling01Dec08

Aug_20_2008 http://www.comlab.ox.ac.uk/projects/chebfun/publications/chebfun_paper.pdf On page 12, it looks like: >> chebpoly(x.^3) Should really be: >> chebpoly(chebfun('x.^3')) ans = 0.250000000000000 0 0.750000000000000 0 Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#Chebpoly20Aug08

Aug_09_2008 Strang and Edelman on Pascal's triangle in matrices. http://www-math.mit.edu/~gs/papers/pascal-work.pdf They report the mapping from x to (x+1), I believe this is the same relationship I use for base conversion. Also, here's a simplified pt() function for octave (and probably matlab). function [x] = pt(Z) % % File pt.m % This function returns Pascal's Triangle in a square upper % triangular matrix with its root at the lower right. % % It can be run as pt(Z), where Z is the size of the square. % % The matrix returned is suitable for polynomial base conversion. % % Transformation matrix, T, used for reordering rows and columns. T = eye(Z); P = abs(pascal(Z,1)); % swap columns of T symmetric around the center. for m = 1:Z/2 TMP=T(:,m); T(:,m) = T(:,Z-m+1); T(:,Z-m+1) = TMP; end % T * P reverses the order of the rows. % TP * T reverses the order of the columns. [x] = T * P * T; Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#PascalAgain09Aug08

Aug_02_2008 So it appears that the reason that multiplication by Pascal's triangle preserves the determinant is because the determinant of Pascal's triangle is "1". It seems that with square matrices, the product of the determinants is the determinant of the product. >> A=rand(30,30); >> B=rand(30,30); >> det(A) * det(B) ans = -0.6344 >> det(A * B) ans = -0.6344 I'm sure I must've known that once upon a time (like 20 years ago), but it had been long ago forgotten. Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#Determinant01Aug08

Jul_31_2008 Looks like matrix multiplication by Pascal's triangle, or it's inverse, preserves the determinant.... scp44@tux64-07:~/programming/matlab$ cat detpres.m SIZE=10; MC=0; % Match Count TOT=0; for lcv=1:10000 A=round(rand(SIZE,SIZE) * 10); B=A * pt(SIZE); C=A * pt(SIZE)^-1; dA=det(A); dB=det(B); dC=det(C); if (dA ~= dB) A B sprintf ('A:B: Determinants are not equal') end if (dA ~= dC) A C sprintf ('A:C: Determinants are not equal') end if ((dA == dB) && (dA == dC)) MC=MC+1; end TOT=TOT+1; end sprintf ('Matched: %d --- Total: %d', MC, TOT) >> detpres ans = Matched: 10000 --- Total: 10000 Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#PascalTriangleDeterminant31July08

Jul_30_2008 This might be something to explore... kinda like matrix factorization(?). >> Z=round(rand(3,3)*10) Z = 10 10 1 2 5 4 10 8 9 >> P=pt(3) P = 1 2 1 0 1 1 0 0 1 >> Pi=P^-1 Pi = 1 -2 1 0 1 -1 0 0 1 >> F=Z*Pi F = 10 -10 1 2 1 1 10 -12 11 >> A=F*P A = 10 10 1 2 5 4 10 8 9 >> Z Z = 10 10 1 2 5 4 10 8 9 >> det(A) ans = 316 >> det(Z) ans = 316 >> svd(A) ans = 21.1999 5.9382 2.5102 >> svd(Z) ans = 21.1999 5.9382 2.5102 >> svd(F) ans = 23.0008 6.1641 2.2288 >> det(F) ans = 316 >> F F = 10 -10 1 2 1 1 10 -12 11 >> Z Z = 10 10 1 2 5 4 10 8 9 Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#MorePascalTriangle30July08

Jul_25_2008 http://home.ccil.org/~remlaps/javascript/GCD.html http://home.ccil.org/~remlaps/weblog_2007/mathematics_index.html#Factors2-13Sep07 I wonder if the factor check method I was playing with could be modified to find greatest common factors? Wow. That was 6 months ago and more? Time flies... And who knows when I'll get time to look into it any more? Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#Checklater25July08

Jun_24_2008 Right to left (ignoring signs) 1st column - constant(1) 2nd column - constant difference (1) 3rd column - 2nd difference (difference of differences) constant (1) 4th column - 3rd difference constant (1). 5th column - 4th difference constant (1). 6th column - 5th difference constant (1). 7th column - 6th difference constant (1). 8th column - *not* 7th (or any) difference constant. Starting points are the 8th row of Pascal's triangle (with alternating signs, starting with positive at left). ans = 128 -448 672 -560 280 -84 14 -1 64 -256 432 -400 220 -72 13 -1 32 -144 272 -280 170 -61 12 -1 16 -80 168 -192 129 -51 11 -1 8 -44 102 -129 96 -42 10 -1 4 -24 61 -85 70 -34 9 -1 2 -13 36 -55 50 -27 8 -1 1 -7 21 -35 35 -21 7 -1 >> [672-432 432-272 272-168 168-102 102-61 61-36 36-21] ans = 240 160 104 66 41 25 15 >> [240-160 160-104 104-66 66-41 41-25 25-15] ans = 80 56 38 25 16 10 >> [80-56 56-38 38-25 25-16 16-10] ans = 24 18 13 9 6 >> [24-18 18-13 13-9 9-6] ans = 6 5 4 3 >> [6-5 5-4 4-3] ans = 1 1 1 Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#ThirdDifference24June08

Jun_24_2008 This is a curious looking matrix... Palindrome with flipped signs across the lower-left to top-right diagonal... Some obvious patterns others that may not be patterns. AT&T doesn't have {21,36,61,102,168,168,272,432,672}. >> pt(8) * pascal(8,1) ans = 128 -448 672 -560 280 -84 14 -1 64 -256 432 -400 220 -72 13 -1 32 -144 272 -280 170 -61 12 -1 16 -80 168 -192 129 -51 11 -1 8 -44 102 -129 96 -42 10 -1 4 -24 61 -85 70 -34 9 -1 2 -13 36 -55 50 -27 8 -1 1 -7 21 -35 35 -21 7 -1 Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#Amatrix24June08

Jun_23_2008 How 'bout conversions from base 10 to bases {1, 0, -1}? # # Convert 1,000 from base 10 to base 1 (sort-of) # >> A=pt(4)^9 A = 1 27 243 729 0 1 18 81 0 0 1 9 0 0 0 1 >> [1 0 0 0] * A ans = 1 27 243 729 # # 1 + 27 + 243 + 729 = 1000 - [x] # === # # Convert 1,000 from base 10 to base 0 (whatever that is) # >> A=pt(4)^10 A = 1 30 300 1000 0 1 20 100 0 0 1 10 0 0 0 1 >> [1 0 0 0] * A ans = 1 30 300 1000 # # 1 * 0^3 = 0 # 30 * 0^2 = 0 # 300 * 0^1 = 0 # 1,000 * 0^0 = 1,000 # ---- Nice! # === # # Convert 1,000 from base 10 to base -1 (is there such a thing? I don't # know how that would work.) # >> A=pt(4)^11 # # Base -1 # A = 1 33 363 1331 0 1 22 121 0 0 1 11 0 0 0 1 >> [1 0 0 0] * A ans = 1 33 363 1331 # # 1331 - 363 + 33 - 1 == 1,000 # So useless though it is, Pascal's triangle seems to be able to deal with polynomials where X = {1, 0, -1} too. Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#WeirdBases23June08

Jun_23_2008 This one works on both matlab and octave, but the abs() manipulation of the cells is a kludge. scp44@tux64-08:~$ cat pt.m function [x] = pt(Z) % % File pt.m % This function returns Pascal's Triangle in a square upper % triangular matrix with its root at the lower right. % % It can be run as pt(Z), where Z is the size of the square. % % The matrix returned is suitable for polynomial base conversion. % % Transformation matrix, T, used for reordering rows and columns. T = zeros(Z); P = pascal(Z,1); for m = 1:Z for n = 1:Z % Put ones in upper-right to lower-left diagonal of T. if ( m + n == Z + 1 ) T(m,n) = 1; end % Get rid of negatives P(m,n) = abs(P(m,n)); end end % T * P reverses the order of the rows. % TP * T reverses the order of the columns. [x] = T * P * T; ### ### octave ### $ echo "pt(9)" | octave === greeting snipped === ans = 1 8 28 56 70 56 28 8 1 0 1 7 21 35 35 21 7 1 0 0 1 6 15 20 15 6 1 0 0 0 1 5 10 10 5 1 0 0 0 0 1 4 6 4 1 0 0 0 0 0 1 3 3 1 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 ### ### matlab ### scp44@tux64-08:~$ echo "pt(9)" | matlab -nodisplay === greeting snipped >> ans = 1 8 28 56 70 56 28 8 1 0 1 7 21 35 35 21 7 1 0 0 1 6 15 20 15 6 1 0 0 0 1 5 10 10 5 1 0 0 0 0 1 4 6 4 1 0 0 0 0 0 1 3 3 1 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#PTkludge23June08

Jun_23_2008 Got it. There's probably a better way, but at least I've got one that works now.... ===== scp44@tux64-07:~$ cat pt.m function [x] = pt(Z) % % File pt.m % This function returns Pascal's Triangle in a square upper % triangular matrix with its root at the lower right. % % It can be run as pt(Z), where Z is the size of the square. % % The matrix returned is suitable for polynomial base conversion. % % Transformation matrix, T, used for reordering rows and columns. T = zeros(Z); % Put one's into the top-right to lower-left diagonal of T. for m = 1:Z for n = 1:Z if ( m + n == Z + 1 ) T(m,n) = 1; end end end % T * P reverses the order of the rows. % TP * T reverses the order of the columns. P = chol(pascal(Z,0),'lower'); [x] = T * P * T; ===== scp44@tux64-07:~$ matlab -nodisplay < M A T L A B > Copyright 1984-2007 The MathWorks, Inc. Version 7.4.0.287 (R2007a) January 29, 2007 To get started, type one of these: helpwin, helpdesk, or demo. For product information, visit www.mathworks.com. >> help pt File pt.m This function returns Pascal's Triangle in a square upper triangular matrix with its root at the lower right. It can be run as pt(Z), where Z is the size of the square. The matrix returned is suitable for polynomial base conversion. >> A=pt(3) A = 1 2 1 0 1 1 0 0 1 >> [1 0 0] * A ans = 1 2 1 >> % Conveniently, 100_10 is 121_9. : -) >> exit ===== scp44@tux64-07:~$ echo "pt(8)" | matlab -nodisplay < M A T L A B > Copyright 1984-2007 The MathWorks, Inc. Version 7.4.0.287 (R2007a) January 29, 2007 To get started, type one of these: helpwin, helpdesk, or demo. For product information, visit www.mathworks.com. >> ans = 1 7 21 35 35 21 7 1 0 1 6 15 20 15 6 1 0 0 1 5 10 10 5 1 0 0 0 1 4 6 4 1 0 0 0 0 1 3 3 1 0 0 0 0 0 1 2 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 >> Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#MatlabFunction-pt-23June08

Jun_23_2008 Working towards automatically generating pascal's triangle with the orientation I like (in matlab). >> [0 0 1;0 1 0;1 0 0] * chol(pascal(3,0),'lower') * [0 0 1;0 1 0;1 0 0] ans = 1 2 1 0 1 1 0 0 1 Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#PascalTriangleMatlab23June08

Jun_22_2008 So far I like octave too... Here's 10000 to base 16 before normalizing digits. octave:1> a=[1 4 6 4 1 > 0 1 3 3 1 > 0 0 1 2 1 > 0 0 0 1 1 > 0 0 0 0 1] a = 1 4 6 4 1 0 1 3 3 1 0 0 1 2 1 0 0 0 1 1 0 0 0 0 1 octave:2> [1 0 0 0 0] * a^-6 ans = 1 -24 216 -864 1296 octave:3> 1 * 16^4 - 24 * 16^3 + 216*16^2 - 864 *16 + 1296 ans = 10000 Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#Octave22June08

Jun_22_2008 I'm liking matlab... I'm also downloading cygwin with octave (supposedly a free matlab clone). Hopefully octave will be this easy too. Here's a polynomial transformation from being in terms of (X) to being in terms of (X+6). A first step of base 10 to base 16 conversion... >> a=[1 3 3 1 0 1 2 1 0 0 1 1 0 0 0 1] a = 1 3 3 1 0 1 2 1 0 0 1 1 0 0 0 1 >> [1 0 0 0] * a^-6 ans = 1 -18 108 -216 >> 1 * 16^3 - 18 * 16^2 + 108 * 16 - 216 ans = 1000 Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#MatlabIsNice22June08

Jun_21_2008 Interesting matrix(?) 1 5 10 10 5 1 1 1 4 6 4 1 1 2 1 3 3 1 1 3 3 1 2 1 1 4 6 4 1 1 1 5 10 10 5 1 Determinant = 0 Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#PascalTriangleAgain21June80

Apr_12_2008 Number of digits needed to represent the Fibonacci numbers in the Fibonacci bases... http://www.postimage.org/image.php?v=Pq1RIvmr I'm imagining a three dimensional matrix where: X-axis is the number (in base 10) Y-Axis is the base (in base 10) Z-Axis is the digit (in base 10) Much of the cube would be empty, but the part that contains information would contain all representations of all Integers. Given a cube like that, I could get to an arbitrary digit of N_b simply by accessing the coordinates. Position (13,2,3) would be 0. The 3rd digit of 13 in base 2. (1101_2) Position (18732,10,4) would be 3. The 4th digit of 18732 in base 10. Position (15,16,1) would be 15. The 1st digit of 15 in base 16 (F_16). Is there a formula to efficiently calculate P(X,Y,Z)? Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#FibReps12Apr08

Apr_09_2008 Number of digits needed to represent a number in a base... Same triangle as before, but easier to view in a spread sheet. http://spreadsheets.google.com/pub?key=pk6AtB0QfKvqluT7ZL5Tydg Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#DigitsNeededTriangle09Apr08

Apr_09_2008 Playing a fantasy investing strategy, I notice that there's another sort of "compound return". If I 1.) invest a principle, P, in some investment, I_0. 2.) wait for it to rise by some percentage, r. 3.) sell P dollars worth of I_0 holding. I'm left with P dollars in cash and r dollars invested in I_0 (minus fees and taxes). The risk of losing money from my remaining investment in I_0 becomes 0. I've already broken even, so that's my worst case scenario. Anything that I_0 does (except going broke) will be profit when I cash it out. I can now proceed back to step 1 with a new investement, I_1. In effect, the same P dollars go to work for me multiple times. Over time, I can wind up with risk-free holdings in {I_0, I_1, I_2, ... I_n}. Of course, the problem is with step 2. There's no guarantee that the investement will ever rise by r percent. I might lose my principle or the investment might do nothing and keep my P dollars tied up indefinitely. I'm still intrigued by this notion though. I think the risk from step 2 can be mitigated by making a number of relatively small sized initial investments. Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#AnotherCompoundInterest09Apr08

Mar_16_2008 http://www.research.att.com/~njas/sequences/A122587 > "Digits 1, 2 and 3 appear alternately and each time in runs whose > lengths are the powers of 4." > > FORMULA: a(n)=floor(n/(4^floor(log[4](n)))) So in general, it looks like the first digit of N in base B will be: a(N,B)=floor(N/(B^floor(log[B](n)))) And here's an implementation of it in emacs lisp... (defun genlog (X B) (/ (log X) (log B))) (defun fd (N B) (floor (/ N (expt B (floor (genlog N B)))))) (fd 39 8) 4 (fd 1321 11) 10 (fd 978432 10) 9 Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#FirstDigitBaseFour16Mar08

Mar_16_2008 http://www.research.att.com/~njas/sequences/A122586 > "Digits 1 and 2 appear alternately and each time in runs whose lengths > are the powers of 3." > > FORMULA --- a(n)=floor(n/(3^floor(log[3](n)))) As seen in the 2nd column of the table below, that series contains the first digits of the base 3 numbers. First Digit of Representation **** Base **** 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 N: ------------------------------------------------------------------------------------ 2: 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3: 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4: 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5: 1, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6: 1, 2, 1, 1, 1, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7: 1, 2, 1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8: 1, 2, 2, 1, 1, 1, 1, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9: 1, 1, 2, 1, 1, 1, 1, 1, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10: 1, 1, 2, 2, 1, 1, 1, 1, 1, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 11: 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 12: 1, 1, 3, 2, 2, 1, 1, 1, 1, 1, 1, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 13: 1, 1, 3, 2, 2, 1, 1, 1, 1, 1, 1, 1, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 14: 1, 1, 3, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 15: 1, 1, 3, 3, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 16: 1, 1, 1, 3, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 17: 1, 1, 1, 3, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 18: 1, 2, 1, 3, 3, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 19: 1, 2, 1, 3, 3, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 20: 1, 2, 1, 4, 3, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 20, 20, 20, 20, 20, 20, 20, 20, 20, 21: 1, 2, 1, 4, 3, 3, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 21, 21, 21, 21, 21, 21, 21, 21, 22: 1, 2, 1, 4, 3, 3, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 22, 22, 22, 22, 22, 22, 22, 23: 1, 2, 1, 4, 3, 3, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 23, 23, 23, 23, 23, 23, 24: 1, 2, 1, 4, 4, 3, 3, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 24, 24, 24, 24, 24, 25: 1, 2, 1, 1, 4, 3, 3, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 25, 25, 25, 25, 26: 1, 2, 1, 1, 4, 3, 3, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 26, 26, 26, 27: 1, 1, 1, 1, 4, 3, 3, 3, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 27, 27, 28: 1, 1, 1, 1, 4, 4, 3, 3, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 28, 29: 1, 1, 1, 1, 4, 4, 3, 3, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#Base3FirstDigits16Mar08

Mar_16_2008 Updated the first-digit square program. Nothing real significant. Simplification and type changes to try to squeeze a little more out of it. Also, I was having small-digit errors from using "float" instead of "double". The updated code is still here: http://home.ccil.org/~remlaps/src/square_benford.cpp.txt I still don't know if the ones/twos ratio is converging or not. It doesn't seem to be, but that's not aesthetically pleasant, so I don't believe it. Maybe I've still got another programming bug or maybe I just haven't been able to get to big-enough numbers to see it converge. I have a real hard time believing that the ratio of [(numbers starting with ones) over (numbers starting with twos)] is going to grow to infinity when ordered like this. That's it for now. Out of time. : -( Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#Sqb16Mar08

Mar_03_2008 Changed the FD array into an int. Now we can see higher numbers. In a 79,999 x 79,999 square, the leading-ones / leading-twos ratio is 2.99616408348. sqb 80000 First Digit Counts: 1: 1605461651 0.762736 2: 535839043 0.254571 3: 268206467 0.127422 4: 161075652 0.0765251 5: 107464688 0.0510552 6: 76822844 0.0364976 7: 57671129 0.0273989 8: 44885448 0.0213245 9: 35936180 0.0170729 10: 29426301 0.0139801 1 / 2: 2.99616408348 2 / 3: 1.99786031246 3 / 4: 1.66509628296 4 / 5: 1.49887049198 5 / 6: 1.39886367321 6 / 7: 1.33208501339 7 / 8: 1.28485131264 8 / 9: 1.24903225899 9 / 10: 1.22122657299 10 / 11: 1.19884729385 Total: 2104872705 Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#AnotherRun03March08

Mar_03_2008 Assuming no bugs (and assuming it's converging), It seems like the one/two ratio is converging very very slowly... == 14000 x 14000 First Digit Counts: 1: 49447222 0.252282 2: 16542519 0.0844006 1 / 2: 2.98909878731 Total: 196000000 == 15000 x 15000 First Digit Counts: 1: 56741853 0.252186 2: 18980887 0.0843595 1 / 2: 2.98942041397 Total: 225000000 == 16000 x 16000 First Digit Counts: 1: 64538619 0.252104 2: 21586200 0.0843211 1 / 2: 2.98980927467 Total: 256000000 Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#FourteenThousandSquare03March08

Mar_02_2008 I found bug in the first digit count program I've been playing with. I was screwing up multiples of the base such as {20, 30, ...}. Here's the revised source code. http://home.ccil.org/~remlaps/src/square_benford.cpp.txt I did these by hand... Table of Representations Number *** Base *** 2 3 4 5 6 7 8 9 A 2 10 2 2 2 2 2 2 2 2 3 11 10 3 3 3 3 3 3 3 4 100 11 10 4 4 4 4 4 4 5 101 12 11 10 5 5 5 5 5 6 110 20 12 11 10 6 6 6 6 7 111 21 13 12 11 10 7 7 7 8 1000 22 20 13 12 11 10 8 8 9 1001 100 21 14 13 12 11 10 9 A 1010 101 22 20 14 13 12 11 10 First *** Count in *** Digit *** Base *** 2 3 4 5 6 7 8 9 10 Total 1 - 9 + 5 + 4 + 5 + 5 + 4 + 3 + 2 + 1 = 38 2 - 0 + 4 + 4 + 2 + 1 + 1 + 1 + 1 + 1 = 15 3 - 0 + 0 + 1 + 1 + 1 + 1 + 1 + 1 + 1 = 7 4 - 0 + 0 + 0 + 1 + 1 + 1 + 1 + 1 + 1 = 6 5 - 0 + 0 + 0 + 0 + 1 + 1 + 1 + 1 + 1 = 5 6 - 0 + 0 + 0 + 0 + 0 + 1 + 1 + 1 + 1 = 4 7 - 0 + 0 + 0 + 0 + 0 + 0 + 1 + 1 + 1 = 3 8 - 0 + 0 + 0 + 0 + 0 + 0 + 0 + 1 + 1 = 2 9 - 0 + 0 + 0 + 0 + 0 + 0 + 0 + 0 + 1 = 1 For comparison with these, which were done programmatically C:\Personal\PROGRA~1\cpp>sqb 10 First Digit of Representation *** Base *** 2 3 4 5 6 7 8 9 10 N: -------------------------- 2: 1, 2, 2, 2, 2, 2, 2, 2, 2, 3: 1, 1, 3, 3, 3, 3, 3, 3, 3, 4: 1, 1, 1, 4, 4, 4, 4, 4, 4, 5: 1, 1, 1, 1, 5, 5, 5, 5, 5, 6: 1, 2, 1, 1, 1, 6, 6, 6, 6, 7: 1, 2, 1, 1, 1, 1, 7, 7, 7, 8: 1, 2, 2, 1, 1, 1, 1, 8, 8, 9: 1, 1, 2, 1, 1, 1, 1, 1, 9, 10: 1, 1, 2, 2, 1, 1, 1, 1, 1, First Digit Counts: 1: 38 0.469136 2: 15 0.185185 3: 7 0.0864198 4: 6 0.0740741 5: 5 0.0617284 6: 4 0.0493827 7: 3 0.037037 8: 2 0.0246914 9: 1 0.0123457 10: 0 0 1 / 2: 2.53333330154 2 / 3: 2.14285707474 3 / 4: 1.16666662693 4 / 5: 1.20000004768 5 / 6: 1.25 6 / 7: 1.33333337307 7 / 8: 1.5 8 / 9: 2 9 / 10: Inf Total: 81 And (assuming no other bugs) here are the first 10 leading-digit ratios in a 2,000 by 2,000 square. First Digit Counts: 1: 1029060 0.257265 2: 347328 0.086832 3: 175186 0.0437965 4: 106015 0.0265037 5: 71537 0.0178843 6: 51420 0.012855 7: 39061 0.00976525 8: 30803 0.00770075 9: 24944 0.006236 10: 20631 0.00515775 First Digit ratios 1 / 2: 2.96279025078 2 / 3: 1.98262417316 3 / 4: 1.65246427059 4 / 5: 1.48196041584 5 / 6: 1.39122915268 6 / 7: 1.31640255451 7 / 8: 1.26809072495 8 / 9: 1.23488616943 9 / 10: 1.20905435085 10 / 11: 1.18623507023 Unsurprisingly, the programming error doesn't seem to have effected the ratios too dramatically. The ratio of leading-ones over leading-twos still seems to be converging just under 3. I also reduced the program run-time from O(N^3 [+ N^2]) to O(N^2) by counting first digits at the same time I calculate them. Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#FirstDigitSquare02March08

Feb_20_2008 First-Digit-1 count / First-digit-2 count, zoomed in. http://www.postimage.org/image.php?v=Pq1BpqU0 I think this is not a well-behaved function. Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#ZoomIn20Feb08

Feb_20_2008 The 5,000 by 5,000 square from 2 through 5,0001 finished. Here are the ratios by first digit. The number of representations of 2 through 5,001 in bases 2 through 5,001 which start with 1 is 2.99383473396 times the number of representations which start with 2. 1 / 2: 2.99383473396 2 / 3: 1.98900818825 3 / 4: 1.65859162807 4 / 5: 1.49368715286 5 / 6: 1.39234983921 6 / 7: 1.32679462433 7 / 8: 1.27852058411 8 / 9: 1.2438287735 9 / 10: 1.21321880817 10 / 11: 1.19107115269 Total: 25000000 I note that the 1/2 ratio is near(but not) 3, the 2/3 ratio is near 2, 3/4 near 1 2/3, 4/5 near 1 1/2 (skip 5/6), 6/7 is near 1 1/3 (skip 7/8) and 8/9 is near 1.25. For whatever that's all worth. Also, I'm thinking that in general, the ratio of 1's and 2's might be calculated as follows: In base 2, all numbers start with 1. In base 3, numbers start with 1 or 2, so (1/2) start with 1 and (1/2) start with 2. In base 4, numbers start with 1, 2 or 3, so (1/3) start with 1 and (1/3) start with 2. etc... So the ratio of first-ones / first-twos in all bases would be: 1 + 1/2 + 1/3 + 1/4 + ... --------------------------- 1/2 + 1/3 + 1/4 + 1/5 + ... And the ratio of first-twos / first-threes in all bases would be: 1/2 + 1/3 + 1/4 + 1/5 + ... --------------------------- 1/3 + 1/4 + 1/5 + 1/6 + ... Unfortunately, so far it's not working that way. The actual tallies don't appear to be converging to the same values as those ratios. It was just an impulse and I've only checked with the windows calculator, so obviously someting needs more thought or less typos. Gotta go. It'll have to wait for later. Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#Five001Reps20Feb08

Feb_20_2008 Here are the ratios of the first 10 adjacent digits in a 4,999 by 4,999 first-digit representation square. 1 / 2: 2.99383974075 2 / 3: 1.98900389671 3 / 4: 1.65859258175 4 / 5: 1.49368917942 5 / 6: 1.39234745502 6 / 7: 1.32678496838 7 / 8: 1.2785333395 8 / 9: 1.24381446838 9 / 10: 1.21322894096 10 / 11: 1.19106757641 Total: 24990001 i.e. Looking at representations of 2 through 5,000 in bases from 2 through 5,000, there are 2.99383974075 times as many representations starting with one as there are starting with two. (assuming no code bugs) Checking 2 through 5,001 now. Probably won't have results until tonight. Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#FiveThousandRatios20Feb08

Feb_20_2008 Here are the ratios of the first 10 adjacent first digit pairs in a 2,000 by 2,000 square.... 1 / 2: 2.99320983887 2 / 3: 1.98440539837 3 / 4: 1.65202283859 4 / 5: 1.4847921133 5 / 6: 1.38785648346 6 / 7: 1.3213917017 7 / 8: 1.26800167561 8 / 9: 1.23386442661 9 / 10: 1.20742189884 10 / 11: 1.18367815018 Got to 2000 by increasing my stack size, as follows... unsigned _stklen = 1048576 * 16; Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#TenRatios20Feb08

Feb_19_2008 Corrected graph... http://www.postimage.org/image.php?v=aV1_RJWi Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#FirstOneTwoRatio19Feb08

Feb_19_2008 Here's a graph I made with the ratio of first-digit one's to first-digit zero's as N gets bigger. http://www.postimage.org/image.php?v=aV1ZbiwS Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#OnesOverTwos19Feb08

Feb_19_2008 The number of representations that begin with 1, in an n-square box is: n Ones 1 1 2 3 3 6 4 9 5 14 6 19 7 25 ... ... AT&T Doesn't have the sequence {1,3,6,9,14,19,25,...} # # Numbers came from the 2nd column of output # $ ./sqb 3 | egrep '1:' 1: 1 1 $ ./sqb 4 | egrep '1:' 1: 3 0.75 $ ./sqb 5 | egrep '1:' 1: 6 0.666667 $ ./sqb 6 | egrep '1:' 1: 9 0.5625 $ ./sqb 7 | egrep '1:' 1: 14 0.56 $ ./sqb 8 | egrep '1:' 1: 19 0.527778 $ ./sqb 9 | egrep '1:' 1: 25 0.510204 # # Note how the size of the square is 2 less than the argument. # the program starts at 2 and runs until 1-less-than the # argument, so 1 and MAX are not included. # $ ./sqb.sparc 3 # 1x1 square 1, First Digit Counts: 1: 1 1 2: 0 0 One / Two: 0 Total: 1 $ ./sqb 4 # 2x2 square 1, 2, 1, 1, First Digit Counts: 1: 3 0.75 2: 1 0.25 3: 0 0 One / Two: Infinity Total: 4 # # Legends are ommitted in the tables above. Here's how # to interpret them... # # First digit # b2 b3 b4... 2: 1, 2, 2, 3: 1, 1, 3, 4: 1, 1, 1, ... ... ... Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#OneReps19Feb08

Feb_18_2008 So the over-all percentage of first-digit 1's seems to be decreasing as the number of representations and size of the base gets bigger. However, the ratio of First-Digit-Ones/First-Digit-Twos seems to be hovering just under 2. My computer won't go much bigger than the examples below with this simple code... The squares are "eye pleasing" too. I looked for 2,2,3,2,3,4,2,3,4,5,2,3,4,5,6,2,3,4,5,6,7,... at AT&T and found this... http://www.research.att.com/~njas/sequences/A131818 It's not quite the same as the upper right hand triangle of the squares below. // // Numbers & bases from 2 through 9. // Number on Y axis, base on X. // $ ./square_benford.x86 10 1, 2, 2, 2, 2, 2, 2, 2, 1, 1, 3, 3, 3, 3, 3, 3, 1, 1, 1, 4, 4, 4, 4, 4, 2, 1, 1, 1, 5, 5, 5, 5, 1, 1, 1, 1, 1, 6, 6, 6, 1, 2, 1, 1, 1, 1, 7, 7, 1, 2, 1, 1, 1, 1, 1, 8, 2, 1, 2, 1, 1, 1, 1, 1, First Digit Counts: 1: 31 0.484375 2: 12 0.1875 3: 6 0.09375 4: 5 0.078125 5: 4 0.0625 6: 3 0.046875 7: 2 0.03125 8: 1 0.015625 9: 0 0 One / Two: 2 Total: 64 // 64 total representations. // 2 times as many 1's as 2's in the first digit. // 48.4375% (31) of all representations start with 1. // // // // Numbers & bases from 2 through 19. // Number on Y axis, base on X. // $ ./square_benford.x86 20 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 2, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 1, 1, 1, 1, 1, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 1, 2, 1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 1, 2, 1, 1, 1, 1, 1, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 2, 1, 2, 1, 1, 1, 1, 1, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 1, 3, 2, 1, 1, 1, 1, 1, 1, 10, 10, 10, 10, 10, 10, 10, 10, 10, 2, 3, 2, 2, 1, 1, 1, 1, 1, 1, 11, 11, 11, 11, 11, 11, 11, 11, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 12, 12, 12, 12, 12, 12, 12, 1, 1, 3, 2, 2, 1, 1, 1, 1, 1, 1, 1, 13, 13, 13, 13, 13, 13, 1, 1, 3, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 14, 14, 14, 14, 14, 1, 1, 3, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 15, 15, 15, 15, 1, 1, 1, 3, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 16, 16, 16, 2, 1, 4, 3, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 17, 17, 1, 1, 4, 3, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 18, 2, 2, 4, 3, 3, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, First Digit Counts: 1: 129 0.398148 2: 46 0.141975 3: 26 0.0802469 4: 18 0.0555556 5: 14 0.0432099 6: 13 0.0401235 7: 12 0.037037 8: 11 0.0339506 9: 10 0.0308642 10: 9 0.0277778 11: 8 0.0246914 12: 7 0.0216049 13: 6 0.0185185 14: 5 0.0154321 15: 4 0.0123457 16: 3 0.00925926 17: 2 0.00617284 18: 1 0.00308642 19: 0 0 One / Two: 1.76923 Total: 324 // 324 total representations. // 1.769 times as many 1's as 2's in the first digit. // 39.8148% (129) of all representations start with 1. // // Numbers & bases from 2 through 1449. // Number on Y axis, base on X. // $ ./square_benford.x86 1450 | pg -n First Digit Counts: 1: 546297 0.26055 2: 182461 0.0870228 3: 92249 0.0439971 4: 56003 0.02671 5: 37733 0.0179963 6: 27313 0.0130266 7: 20722 0.00988313 8: 16390 0.00781703 ... 1439: 10 4.76939e-06 1440: 9 4.29245e-06 1441: 8 3.81551e-06 1442: 7 3.33857e-06 1443: 6 2.86163e-06 1444: 5 2.3847e-06 1445: 4 1.90776e-06 1446: 3 1.43082e-06 1447: 2 9.53878e-07 1448: 1 4.76939e-07 1449: 0 0 One / Two: 1.97792 Total: 2096704 // 2096704 total representations. // 1.97792 times as many 1's as 2's in the first digit. // 26.055% (546297) of all representations start with 1. /*********************************\ * And here's the C++ source code * \*********************************/ $ cat square_benford.cpp #include #include using namespace std; int main(int argc, char* argv[]) { int MAX=10; if (argc == 2 ) MAX=atoi(argv[1]); int FD[MAX-2][MAX-2]; // Number for ( int y=2; y x ) DIV = DIV/x; if ( DIV == 0 ) { FD[x-2][y-2] = MOD; if ( MAX < 31 ) cout << MOD << ", "; } else if ( MOD != 0 ) { FD[x-2][y-2] = DIV; if ( MAX < 31 ) cout << DIV << ", "; } else { FD[x-2][y-2] = 1; if ( MAX < 31 ) cout << 1 << ", "; } } if ( MAX < 31 ) cout << "\n"; } cout << "First Digit Counts:\n"; int count[MAX-2], total=0; for (int lcv=1; lcv Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#FirstDigitAllReps18Feb08

Feb_18_2008 Following up on this entry... http://home.ccil.org/~remlaps/weblog_2007/mathematics_index.html#BenfordCont011406 If I did it right in my head, here are all representations for numbers from 2 through 9 in bases 2 through 9. 2 - 10, 2, 2, 2, 2, 2, 2, 2 3 - 11, 10, 3, 3, 3, 3, 3, 3 4 - 100, 11, 10, 4, 4, 4, 4, 4 5 - 101, 12, 11, 10, 5, 5, 5, 5 6 - 110, 20, 12, 11, 10, 6, 6, 6 7 - 101, 21, 13, 12, 11, 10, 7, 7 8 - 1000, 22, 20, 13, 12, 11, 10, 8 9 - 10001, 100, 21, 14, 13, 12, 11, 10 So out of 64 representations, 33 of them start with "1". 1 - 31 2 - 12 3 - 6 4 - 5 5 - 4 6 - 3 7 - 2 8 - 1 9 - 0 So I guess it's really not surprising that so many numbers in random data start with "1". It appears that the a high proportion of the possible representations start with 1 (when counted this way). I think if I find some time, I'll program this and see what happens when the square gets bigger. http://home.ccil.org/~remlaps/weblog_2007/mathematics_index.html#BenfordLaw011207 As I saw 50% of the representations of 4,000 began with 1, I'm expecting the trend will stay more or less the same as the numbers get larger. Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#Benford18Feb08

Jan_26_2008 Here are the numbers from 1 to 100, ordered randomly and split into groups of 10. Grading on a curve, the 10 "low performers" are: {4, 0, 2, 14, 1, 10, 19, 3, 6, 5} and the 10 "high performers" are: {87, 95, 75, 98, 91, 82, 83, 99, 89, 90}. 7 86 73 79 44 16 4 32 49 87 - 4 - 87 40 39 76 0 64 46 88 95 53 15 - 0 - 95 37 69 11 2 78 42 62 21 65 75 - 2 - 75 51 14 92 66 97 58 96 29 98 56 - 14 - 98 91 20 23 70 1 9 84 68 57 55 - 1 - 91 10 18 74 27 82 54 36 26 25 60 - 10 - 82 41 61 19 48 67 83 34 45 43 47 - 19 - 83 93 94 59 50 33 24 3 77 99 31 - 3 - 99 52 30 81 89 6 22 28 38 17 71 - 6 - 89 85 13 12 80 35 72 5 8 90 63 - 5 - 90 So the top 10 picks were 50/50 and the bottom 10 were 70/30. 60% right and 40% wrong for the 20 values at the boundaries. I'm thinking this is what "grading on a curve" does to rankings. Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#TopNum26Jan08

Jan_25_2008 Rumor has it that some companies employ an evaluation policy based upon statistics. It goes something like this. Say there are 5 possible rankings. 1-Lousy 2-Bad 3-Mediocre 4-Good 5-Excellent And say a manager has 20 employees and (s)he thinks they're all excellent. HR comes along and says "no you don't have 20 excellent employees! I don't care what you think. The talent of your employees follows the normal curve. Rerank them. You've got 2 who are excellent, 2 who are lousy, 3 good, 3 bad and 10 average." Now basically, this is just the old "grading on a curve" that we all remember from college and high school. When I was in school, it seemed tough but fair, and sometimes it even worked in my favor. What's really going on though, I wonder. The first thing that jumps out at me is that the Normal curve should apply to a large population. I don't think any statistician would really expect a small sample, like 20, to follow a normal curve. So if our manager works in a company with say 10,000 employees who's skills are distributed according to a Normal curve, there might very well be a group of 20 who are all excellent and another group of 20 who are all lousy. Now along comes HR, instructing the managers of both groups to re-rank. Now we've got 3 excellent employees ranked as lousy and 3 lousy ones ranked as excellent. To say nothing of the other, intermediate ranks. I'm thinking that company-wide, this method of ranking is going to get almost as many evaluations wrong as it gets right. Maybe even more. The second thing that jumps out at me is what happens at the management and executive levels. If we start with the assumption that people are managers for a reason and therefore, they are probably good or excellent in reality, then the HR policy is certain to mis-rank the majority of managers. If I ever get time, it might be interesting to figure out exactly how bad the expected mis-ranking would be under a system like that. I am now thinking that grading on a curve, for teachers, is probably a really bad idea. Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#NormalCurveHarmful25Jan08

Jan_20_2008 Base case: I can find 1 second every day to ride the exercise bike. Inductive hypothesis: If I can find k seconds to ride the exercise bike every day, then I can find (k+1) seconds to ride the exercise bike every day. Discussion: Somewhere between {1 second} and {23 hours, 59 minutes and 59 seconds}, that argument goes from being true to being false. I am pretty sure I could find one second every day to ride the exercise bike. I'm pretty sure I can't find 24 hours every day for exercise. Where does the argument break down? Which second makes it false? Maybe the cost/weight of an additional second grows exponentially and maybe recursive proofs aren't appropriate for finite sets or exponential growth. Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#RecursiveProof20Jan08

Jan_06_2008 I don't remember if I did this before or not. I don't think so. I've been sidetracked for a while. Here's a lisp implementation of the general mod routine I mentioned a few months back. The mod algorithm was mentioned here- http://home.ccil.org/~remlaps/weblog_2007/mathematics_index.html#Factors2-13Sep07 So we can check whether a base 10 number (F_10) is a factor of a number in an arbitrary base (N_b), or "what is mod(N_b, F_10)?" without converting either to a common base. (defun gmod (lastRem digits base factor) "*** *** Take (digits) in ~base~ and get mod ~factor~ (base 10) *** using base 10 arithmetic. No need to convert digits *** or factor to a common base. ***" (cond ((equal (first digits) nil) nil) ;;; Bad place to be. ((equal (rest digits) nil) ;;; One digit (mod (+ lastRem (first digits)) factor)) ;;; Recursive call for multiple digits. (t (gmod ;;; nested lastRem gets filled with this mod result ;;; Consuming previous remainder and next digit. (mod (* (+ lastRem (first digits)) base) factor) (rest digits) base factor)))) ;;; ;;; If gmod gives back 0, then factor is a factor of the number given by (digits) ;;; in base ~base~. ;;; (defun genmod (DList base factor) "Calling function for gmod so I don't have to supply the leading 0" (gmod 0 DList base factor)) (defun isFactor (DList base factor) "Is ~factor~ a factor of DList in base ~base~?" (equal 0 (gmod 0 DList base factor))) And here's a sample run c:\Personal\programming\lisp>gcl1 gcl1 GCL (GNU Common Lisp) (2.6.1) Wed Oct 8 14:58:05 EAST 2003 Licensed under GNU Library General Public License Dedicated to the memory of W. Schelter Use (help) to get some basic information on how to use GCL. >(require 'modfun.lisp) Loading modfun.lisp Finished loading modfun.lisp T >(genmod '(4 91 37) 100 17) ;;; Base 100 7 >(mod 49137 17) ;;; Compare to standard mod. 7 >(genmod '(37107) 100000 31) ;;; Base 100,000 0 >(genmod '(37 107) 1000 31) ;;; Base 1,000 0 >(genmod '(3 71 07) 100 31) ;;; Base 100 0 >(genmod '(3 7 1 0 7) 10 31) ;;; Base 10 0 >(mod 37107 31) ;;; Standard mod for comparison 0 >(genmod '(37 108) 1000 31) ;;; A related example for a non-factor 1 >(genmod '(63 1 23 4) 62 131) ;;; Base 62 2 >(mod (+ (* 63 62 62 62) (* 1 62 62) (* 23 62) 4) 131) ;;; Standard mod comparison 2 Link for this entry: http://home.ccil.org/~remlaps/weblog_2008/mathematics_index.html#Genmod06Jan07