Message Boards

GROUPS:

Staff Picks Mathematics Algebra Discrete Mathematics External Programs and Systems Wolfram Language Symbolic Computations Computer-Based Maths

Shenghui Yang

Using RunProcess to study lunar arithmetic on Numberphile

Shenghui Yang, WOLFRAM

Posted 4 months ago

1455 Views

2 Replies

6 Total Likes

Follow this post

lunarprimes

To understand the background and rule of the lunar arithmetic, please check the following Numberphile video and research paper (by David Applegate)

◼

YouTube link: https://youtu.be/cZkGeR9CWbk

◼

Paper: http://emis.impa.br/EMIS/journals/JIS/VOL14/Sloane/carry2.pdf

The algorithms is directly from David Applegate’s C code

◼

Link to the source code: https://oeis.org/A087061/a087061.txt

Build the c code on local machine with gcc (MacOS/Linux/Cygwin for windows)

gcc -o dismal_eis dismal_eis.c

Because this code is implemented with stdin/stdout, I can stream the output directly into WL with RunProcess function. Everything I write in the notebook is nothing more than a simple wrapper to parse and reformat the textual result from the C program. In this way you can save lots of time if you just want to play around the existent code with minimal coding effort, in particular, to implement the divisor finder and prime test function.

Add the location of the runtime

In[]:=

$LunarArithmeticRuntimePath="C:\\cygwin64\\home\\sheng\\dismal_eis.exe";

(for Cygwin you need to copy and paste cygwin1.dll to same directory as the runtime)

Test if the runtime works properly. If nothing wrong with compilation, I can launch the help info with no problem

In[]:=

RunProcess[{$LunarArithmeticRuntimePath}]["StandardError"]

Out[]=

Usage: /sheng/dismal_eis cmd [args...]available cmds are: div a b Outputs c such that b*c = a divs a Outputs b such that b*c = a for some c add a b Outputs a+b mul a b Outputs a*b pcount a b Counts the primes in [a,b] dinfo a Outputs info about divisors of a dinfo a b Outputs info about divisors of n, a<=n<=b pdinfo a Outputs info about prime divisors of a pdinfo a b Outputs info about prime divisors of n, a<=n<=b help_dinfo Describes the output of dinfo and pdinfo help Outputs this text

Build formal operator based on

In[]:=

CircleTimes[a_Integer,b__Integer]:=Fold[LunarTimes,a,{b}]

In[]:=

CirclePlus[a_Integer,b__Integer]:=Fold[LunarPlus,a,{b}]

In[]:=

Table[Row[{i," ⊗ ",i," = ",i⊗i}],{i,1,9}]//Column

Out[]=

1 ⊗ 1 = 1

2 ⊗ 2 = 2

3 ⊗ 3 = 3

4 ⊗ 4 = 4

5 ⊗ 5 = 5

6 ⊗ 6 = 6

7 ⊗ 7 = 7

8 ⊗ 8 = 8

9 ⊗ 9 = 9

Some two-digit square table

In[]:=

Table[Row[{i," ⊗ ",i," = ",i⊗i}],{i,10,50}]//Column

Out[]=

10 ⊗ 10 = 100

11 ⊗ 11 = 111

12 ⊗ 12 = 112

13 ⊗ 13 = 113

14 ⊗ 14 = 114

15 ⊗ 15 = 115

16 ⊗ 16 = 116

17 ⊗ 17 = 117

18 ⊗ 18 = 118

19 ⊗ 19 = 119

20 ⊗ 20 = 200

21 ⊗ 21 = 211

22 ⊗ 22 = 222

23 ⊗ 23 = 223

24 ⊗ 24 = 224

25 ⊗ 25 = 225

26 ⊗ 26 = 226

27 ⊗ 27 = 227

28 ⊗ 28 = 228

29 ⊗ 29 = 229

30 ⊗ 30 = 300

31 ⊗ 31 = 311

32 ⊗ 32 = 322

33 ⊗ 33 = 333

34 ⊗ 34 = 334

35 ⊗ 35 = 335

36 ⊗ 36 = 336

37 ⊗ 37 = 337

38 ⊗ 38 = 338

39 ⊗ 39 = 339

40 ⊗ 40 = 400

41 ⊗ 41 = 411

42 ⊗ 42 = 422

43 ⊗ 43 = 433

44 ⊗ 44 = 444

45 ⊗ 45 = 445

46 ⊗ 46 = 446

47 ⊗ 47 = 447

48 ⊗ 48 = 448

49 ⊗ 49 = 449

50 ⊗ 50 = 500

I use Fold function in the definition of the operators, therefore both lunar plus and lunar multiplication supports more than 3 operands

In[]:=

48⊗49⊗51⊗52

Out[]=

44521

Associative, Commutative and Distributive

Lunar arithmetic is basically mixture of max and min operation. It is suffice to show if the following examples holds, all properties holds for lunar arithmetic

◼

Associative property

In[]:=

Max[a,Max[b,c]]===Max[Max[a,b],c]

Out[]=

True

In[]:=

Min[a,Min[b,c]]===Min[Min[a,b],c]

Out[]=

True

◼

Commutative property

In[]:=

Max[a,b]Max[b,a]

Out[]=

True

In[]:=

Min[a,b]Min[b,a]

Out[]=

True

◼

Distributive property

In[]:=

Resolve[ForAll[{a,b,c},Min[a,Max[b,c]]Max[Min[a,b],Min[a,c]]],Reals]

Out[]=

True

No subtraction. No division.

Because result from the inverse operation of lunar plus and lunar division are not unique. For instance

In[]:=

8⊗8

Out[]=

In[]:=

9⊗8

Out[]=

So if I define lunar division, 8 lunar dividing 8 returns either 9 or 8. Same argument applies to the inverse operation of lunar plus. But we can still call 9 and 8 divisors of 8.

Multiplicative Identity

On the moon, 9 is one

In[]:=

9⊗1234

Out[]=

1234

Same as

In[]:=

Min[9,#]&/@IntegerDigits[1234]

Out[]=

{1,2,3,4}

In[]:=

FromDigits[%]

Out[]=

1234

Some number theory functions

Display help information for divisors of a

In[]:=

RunProcess[{$LunarArithmeticRuntimePath,"help_dinfo"}]["StandardOutput"]

Out[]=

dinfo and pdinfo output a header row, and then a row for eachnumber requested. Each row consists of '|' separated fields.The fields are: n the number 9ish 1 if n contains a 9, 0 otherwise prime 1 if n is prime, 0 otherwise n is prime if n != 9, and if n=a*b, then either a==9 or b==9 pseudoprime 1 if n has a divisor with length in [2,n) divisors the number of divisors of n bd_divisors the number of divisors of n in [1,n] dsum_divisors the dismal sum of the divisors of n dsum_bd_divisors the dismal sum of the divisors of n in [1,n] dsum_bd_divisors2 the dismal sum of the divisors of n in [1,n) dsum_ne_divisors the dismal sum of the divisors of n that are != n 2*n the dismal product 2*n sum_divisors the normal sum of the divisors of n sum_bd_divisors the normal sum of the divisors of n in [1,n] sum_bd_divisors2 the normal sum of the divisors of n in [1,n) sum_ne_divisors the normal sum of the divisors of n that are != nIn addition, pdinfo includes the following fields about prime divisors prime_divisors the number of prime divisors of n sum_prime_divisors the dismal sum of the prime divisors of n prod_prime_divisors the dismal product of the prime divisors of n because numbers like 111...111 have very many prime divisors, this product is capped at length 16*length(n)

Test some number from https://oeis.org/A087097

In[]:=

LunarPrimeQ[109]

Out[]=

True

Check the information of the divisors of a given number. Both 19 and 419 are primes and the third entry displays one.

In[]:=

LunarPrimeDivsorsInfo[19]//Dataset

Out[]=

n	19
9ish	1
prime	1
pseudoprime	1
divisors	2
bd_divisors	2
dsum_divisors	19
dsum_bd_divisors	19
dsum_bd_divisors2	9
dsum_ne_divisors	9
2*n	12
sum_divisors	28
sum_bd_divisors	28
sum_bd_divisors2	9
sum_ne_divisors	9
prime_divisors	1
sum_prime_divisors	19
prod_prime_divisors	19

In[]:=

LunarPrimeDivsorsInfo[419]//Dataset

Out[]=

n	419
9ish	1
prime	1
pseudoprime	1
divisors	2
bd_divisors	2
dsum_divisors	419
dsum_bd_divisors	419
dsum_bd_divisors2	9
dsum_ne_divisors	9
2*n	212
sum_divisors	428
sum_bd_divisors	428
sum_bd_divisors2	9
sum_ne_divisors	9
prime_divisors	1
sum_prime_divisors	419
prod_prime_divisors	419

10 is not a prime. So I have the third entry says prime = 0

In[]:=

LunarPrimeDivsorsInfo[10]//Dataset

Out[]=

n	10
9ish	0
prime	0
pseudoprime	1
divisors	18
bd_divisors	10
dsum_divisors	99
dsum_bd_divisors	19
dsum_bd_divisors2	9
dsum_ne_divisors	99
2*n	10
sum_divisors	495
sum_bd_divisors	55
sum_bd_divisors2	45
sum_ne_divisors	485
prime_divisors	1
sum_prime_divisors	90
prod_prime_divisors	90

Find smallest prime. No lunar prime available from 1 to 10

In[]:=

Cases[Range[10],_?LunarPrimeQ]

Out[]=

{}

for example 8 is not prime because

In[]:=

8⊗8

Out[]=

the first one occurs is 19

In[]:=

Cases[Range[20],_?LunarPrimeQ]

Out[]=

{19}

In[]:=

DiscretePlot[LunarPrimePi[n],{n,1,1000}]

Out[]=

Same as

In[]:=

ListPlot[Accumulate[Boole[LunarPrimeQ/@Range[1000]]]]

Out[]=

Try

In[]:=

LunarPrimeQ[1000]

Out[]=

False

In[]:=

LunarPrimePi[1000]

Out[]=

Lunar Prime’s Rarity

Compare the number of lunar primes and regular primes below 100.

In[]:=

lp=Cases[Range[100],_?LunarPrimeQ]

Out[]=

{19,29,39,49,59,69,79,89,90,91,92,93,94,95,96,97,98,99}

In[]:=

p=Cases[Range[100],_?PrimeQ]

Out[]=

{2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97}

Visualize on number line

NumberLinePlot[{lp,p}]

Out[]=

Find the numbers are both prime and lunar prime under 1000

In[]:=

With[{list=Range[1000]},Intersection[Cases[list,_?LunarPrimeQ],Cases[list,_?PrimeQ]]]

Out[]=

{19,29,59,79,89,97,109,409,419,439,509,619,659,709,719,739,769,809,829,839,859,907,919,929,937,947,967}

In[]:=

Length[%]/1000.

Out[]=

0.002

In[]:=

Once[pl2=(LunarPrimePi[#]/N[#])&/@(10^Range[2,7])]

Out[]=

{0.18,0.099,0.1638,0.22095,0.264312,0.315911}

In[]:=

p2=PrimePi[(10^Range[2,7])]/(10.^Range[2,7])

Out[]=

{0.25,0.168,0.1229,0.09592,0.078498,0.0664579}

As my range of test increases, there are much more lunar primes than primes. This is because

In[]:=

ListLogLinearPlot[{Transpose[{10^Range[2,7],pl2}],Transpose[{10^Range[2,7],p2}]},JoinedTrue,PlotLegends{"Lunar prime density","prime density"}]

Out[]=

	Lunar prime density
	prime density

Unlimited Supply of Lunar Primes

In the last section of the YouTube video Sloan constructs the following sequence of numbers

In[]:=

LunarPrimeQ/@{19,109,1009,10009,100009}

Out[]=

{True,True,True,True,True}

Because such list in unlimited, the number of lunar primes are infinity.

Function Definition

In[]:=

$LunarArithmeticRuntimePath="C:\\cygwin64\\home\\sheng\\dismal_eis.exe";

In[]:=

LunarTimes[a_Integer,b_Integer]:=ToExpression@RunProcess[{$LunarArithmeticRuntimePath,"mul",ToString[a],ToString[b]}]["StandardOutput"]

In[]:=

LunarPlus[a_Integer,b_Integer]:=ToExpression@RunProcess[{$LunarArithmeticRuntimePath,"add",ToString[a],ToString[b]}]["StandardOutput"]

In[]:=

LunarPrimeDivsorsInfo[a_Integer]:=AssociationThread@@Fold[StringSplit,RunProcess[{$LunarArithmeticRuntimePath,"pdinfo",ToString[a]}]["StandardOutput"],{"\n","|"}]

In[]:=

LunarPrimePiInfo[n_Integer]:=RunProcess[{$LunarArithmeticRuntimePath,"pcount","1",ToString[n]}]["StandardOutput"]

In[]:=

LunarPrimePi[n_Integer]:=ToExpression[Fold[StringSplit,LunarPrimePiInfo[n],{"\n",WhitespaceCharacter}][[-1,1]]]

In[]:=

LunarPrimeInfo[n_Integer]:=AssociationThread@@Fold[StringSplit,RunProcess[{$LunarArithmeticRuntimePath,"pdinfo",ToString[n]}]["StandardOutput"],{"\n","|"}]

In[]:=

LunarPrimeQ[n_Integer]:=AssociationThread[Rule@@Fold[StringSplit,LunarPrimeInfo[n],{"\n","|"}]]["prime"]/.{"1"True,"0"False}

POSTED BY: Shenghui Yang

Answer

2 Replies

Sort By:

Moderation Team

Moderation Team, WOLFRAM

Posted 4 months ago

-- you have earned *Featured Contributor Badge* Your exceptional post has been selected for our editorial column *Staff Picks* http://wolfr.am/StaffPicks and Your Profile is now distinguished by a *Featured Contributor Badge* and is displayed on the Featured Contributor Board. Thank you!

POSTED BY: Moderation Team

Answer

Moderation Team

Moderation Team, WOLFRAM

Posted 2 months ago

Moderator Note: this post was highlighted on the Wolfram's official social media channels. Thank you for your contribution. We are looking forward to your future posts. Twitter Facebook

POSTED BY: Moderation Team

Answer

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback