OK, and the solution you're comparing it to is when you use primepi(x) to determine exactly how many primes you're going to have, right?

There are ways around that, although they're not as nice. For example, if you have a reasonable upper bound, you can do v=vector(upperbound(x)) and then trim off the remaining values with a technique like I mentioned before. Faiing that, you can resize the vector periodically, then trim the values. If that's too much hassle, you can build up the vector by concatenating, but this is inefficient if the time spent per element is small compared to the number of elements.

[CODE]v=vector(1000,n,if(isprime(n),n,0));c=vector(1,n,0);for(i=1,1000,if(v[i]!=0,c=concat(c,v[i])));for(i=1,(#c)-1,c[i]=c[i+1]);c=vector(#c-1,i,c[i])[/CODE]

this what i came up with on your advice earlier on in the thread.

and your pm.

1/3 times I tested each my 2 idea gave 47 ms for the latter it gave me 46 47 and 63 if you average these 3 you get it's slower I know why to get the same output you need to remove all the 0's so it takes more operations. which take longer especially if my fan starts running. so my first was faster and could store more primes as it didn't waste a entry of a 0 each time it wasn't prime hence 8771 0 removals in the extra time. I tested until 10000 for each.

Your code, slightly generalized:
[code]sm(lim)={
my(v=vector(lim,n,if(isprime(n),n,0)),c=vector(1,n,0));
for(i=1,lim,
if(v[i]!=0,c=concat(c,v[i]))
);
for(i=1,(#c)-1,
c[i]=c[i+1]
);
vector(#c-1,i,c[i])
};[/code]

My code:
[code]crg(lim)={
my(v=vector(primepi(lim)),i=0);
forprime(p=2,lim,v[i++]=p);
v
};[/code]

Your code takes several minutes (hasn't finished by the time I posted this) to give the primes up to a million. Mine takes 187 milliseconds. The built-in primes() command is almost instantaneous, maybe 1 millisecond.

I mean the code I based on what you said versus my original with the 2 change.

on mine the first can be as much as 33% slower.

Okay. So you are busy accounting for the vector portion of the project?

It begins: (My idea):

knfacsieve(a, b, c, d, pmax) = {
for(n=a, b,
forprime(p=2,pmax,...

I just can't finish it up. I think what I posted above is simply going to repeat the forprime loop n times.

I think I have the appropriate idea at hand about what the sieve is supposed to do, but am unable to write the script to make it work.

what is it in basic terms I could look up easily I'll try.

[QUOTE=science_man_88]what is it in basic terms I could look up easily I'll try.
[/QUOTE]

It is a sieve that kicks out bad k-values for numbers of the form k * b![sup]n[/sup] + 1 which are divisible by small primes, where the user is allowed to choose the k-range, the factorial base and the exponent, as well as the largest prime it will sieve up to. The output will be written to a file for testing. The user will have to make the modifications to ensure that the range will be testable, preferably with a decent or advanced text editor that is capable of doing so.

Where k < b![sup]n[/sup], and where k, b, and n are positive integers.

That's pretty much what I plan to do. NewPGen will not accept any base larger than 10 digits, apparently. Hopefully it can be done.

I will then compare this to a script that trial-divides each candidate to 1048576, which is rather easy to construct. (Or 2[sup]20[/sup])

Here's an example for k * 19![sup]6[/sup] + 1:
[code]factor is 241
factor is 579637
factor is 710341
factor is 31249
factor is 67
factor is 169937
factor is 47
factor is 20269
factor is 489449
factor is 43
factor is 643
factor is 3907
factor is 2897
factor is 35993
factor is 103
factor is 613
factor is 733
factor is 23
factor is 919
factor is 2659
factor is 32377
factor is 324529
factor is 3677
factor is 1303
factor is 3583
factor is 18959
factor is 1005427
factor is 164621
factor is 1949
factor is 104479
factor is 31
factor is 101
factor is 29
factor is 41
factor is 14281
factor is 223
factor is 37
factor is 829
factor is 1039
factor is 1030951
factor is 271
factor is 373
factor is 1289
factor is 47041
factor is 318403
factor is 449
factor is 233
factor is 699649
factor is 711793
factor is 23
factor is 5827
factor is 1471
factor is 53
factor is 191
factor is 349
factor is 41177
factor is 157
factor is 947
factor is 61
factor is 50591
factor is 59
factor is 239
factor is 142169
factor is 113
factor is 42443
factor is 172801
factor is 97
factor is 197
factor is 1721
factor is 43
factor is 1217
factor is 47
factor is 21179
factor is 383
factor is 29
factor is 31
factor is 599
factor is 23
factor is 83
factor is 797
factor is 21169
factor is 33829
factor is 1327
factor is 79
factor is 229
factor is 937
factor is 6883
factor is 37
factor is 73
factor is 1663
factor is 41
factor is 863
factor is 6121
factor is 165047
factor is 181
factor is 1277
factor is 167
factor is 751
factor is 26669
factor is 193
factor is 67
factor is 71
factor is 151
factor is 433
factor is 3271
factor is 461
factor is 139
factor is 1579
factor is 137
factor is 335077
factor is 13441
factor is 23
factor is 18307
factor is 89
factor is 577
factor is 109
factor is 199
factor is 7417
factor is 29
factor is 263
factor is 31
factor is 7333
factor is 571
factor is 117937
factor is 311
factor is 53
factor is 1259
factor is 329717
factor is 71663
factor is 63709
factor is 43
factor is 5683
factor is 127
factor is 65957
factor is 311713
factor is 10091
factor is 107
factor is 47
factor is 163
factor is 37
factor is 59
factor is 61
factor is 17497
factor is 224633
factor is 23
factor is 563
factor is 41
factor is 10939
factor is 397
factor is 4153
factor is 2843
factor is 29
factor is 5821
factor is 103
factor is 659
factor is 859
factor is 331
factor is 4079
factor is 873139
factor is 31
factor is 10909
factor is 1571
factor is 207113
factor is 131
factor is 149
factor is 293
factor is 929
factor is 5281
factor is 22973
factor is 101
factor is 421
factor is 23
factor is 194609
factor is 4243
factor is 457
factor is 6659
factor is 89983
factor is 209203
factor is 4951
factor is 4423
factor is 14629
factor is 11131
factor is 43
factor is 67
factor is 73
factor is 37
factor is 1013
factor is 34721
factor is 131711
factor is 21767
factor is 110479
factor is 79
factor is 29
factor is 71
factor is 83
factor is 251
factor is 1613
factor is 14549
factor is 25933
factor is 53
factor is 479
factor is 12347
factor is 52529[/code]

[CODE]sievex(kmin,kmax,b,n,pmax)= for(k=kmin,kmax,forprime(p=1,pmax,if((k*b!^n+1)%p==0,,print("no factor found under "pmax" for "k"*"b"!^"n"+1"))))[/CODE]

best i could do i know it prints a lot I messed something up but it what I could do.

I know what I messed up now it's that I didn't check that all factors were checked for each.

[QUOTE=science_man_88]best i could do i know it prints a lot I messed something up but it what I could do.
[/QUOTE]

You asked it for for every candidate that is not divisible by a certain prime. It is repeating the loop over and over for every prime.

[QUOTE=science_man_88;226740][CODE]sievex(kmin,kmax,b,n,pmax)= for(k=kmin,kmax,forprime(p=1,pmax,if((k*b!^n+1)%p==0,,print("no factor found under "pmax" for "k"*"b"!^"n"+1"))))[/CODE]

best i could do i know it prints a lot I messed something up but it what I could do.[/QUOTE]

[CODE]sievex(kmin,kmax,b,n,pmax)= for(k=kmin,kmax,forprime(p=1,pmax,if((k*b!^n+1)%p==0 && (k*b!^n+1)!=p,print("found factor " p " for "k"*"b"!^"n"+1"))))[/CODE]

give it a wild limit and it prints a lot.

I think I could improve it one more way. should i try to help it more ?

[QUOTE=science_man_88]I think I could improve it one more way. should i try to help it more ?
[/QUOTE]

Do what you can.

would you rather it spit out when it found a factor or when it hasn't(less printing I think).

one thing to stop limitation (though going far enough with pmax it may slow down)
a code to check if it's above primelimit if not use the code already made if not use another code. so basically

[CODE]if(pmax>default(primelimit),code1,code2)[/CODE]

code one could increase primelimit and if it can't go high enough switch to code 2 maybe. might be harder to code.

[CODE]sievex(kmin,kmax,b,n,pmax)={a=0;for(k=kmin,kmax,forprime(p=1,pmax,if((k*b!^n+1)%p==0 && (k*b!^n+1)!=p,a=a+1)));if(a==0,print("no factor found for "k"*"b"!^"n"+1 in this range");a=0)}[/CODE]

can't get this to work like I wanted.

[QUOTE=science_man_88]can't get this to work like I wanted.
[/QUOTE]

What seems to be the error?

[QUOTE=3.14159;226739]Here's an example for k * 19![sup]6[/sup] + 1:
[code]
k=1: factor is 241
k=1: factor is 579637
k=1: factor is 710341
k=2: factor is 31249
k=3: factor is 67
k=4: factor is 169937
k=5: factor is 47
k=7: factor is 20269
k=7: factor is 489449
k=9: factor is 43
k=9: factor is 643
k=9: factor is 3907
k=9: factor is 2897
[/code][/QUOTE]

I've inserted the k-value where the factors come from as example and verification (see [url=http://factordb.com/search.php?query=n*%2819%21%5E6%29%2B1&v=n&n=1&EC=1&E=1&Prp=1&P=1&C=1&FF=1&CF=1&of=H&pp=50&sw=Update]here[/url] for more factorizations. The first PRP is for k=150!).

[QUOTE=kar_bon]I've inserted the k-value where the factors come from as example and verification (see here for more factorizations. The first PRP is for k=150!).
[/QUOTE]

The first PRP = [code]57133839564458545904789328652610540031895535786011264182548375833179829124845398393126574488675311145377107878746854204162666250198684504466355949195922066574942592095735778929325357290444962472405416790722118445437122269675520000000000000000000000000000000000000 * 19![sup]6[/sup] + 1[/code]? O rly?

your guess is as good as mine there Pi:

[CODE](20:53) gp > sievex(1,10,6,1,500000)
found factor 7 for 1*6!^1+1
found factor 103 for 1*6!^1+1
found factor 11 for 2*6!^1+1
found factor 131 for 2*6!^1+1
found factor 43 for 4*6!^1+1
found factor 67 for 4*6!^1+1
found factor 13 for 5*6!^1+1
found factor 277 for 5*6!^1+1
found factor 29 for 6*6!^1+1
found factor 149 for 6*6!^1+1
found factor 71 for 7*6!^1+1
found factor 7 for 8*6!^1+1
found factor 823 for 8*6!^1+1
found factor 19 for 10*6!^1+1
found factor 379 for 10*6!^1+1
time = 281 ms.
<)%p==0 && (k*b!^n+1)!=p,a=a+1)));if
time = 0 ms.
%270 = (kmin,kmax,b,n,pmax)->a=0;for
(20:54) gp > sievex(1,10,6,1,500000)[/CODE]

first sieve gives the result shown first the second sieve should show no factor for 3*6!^1+1 as it only prints if it finds no factor but it hasn't printed anything. both test had the same sievex() values

[QUOTE=3.14159;226752]The first PRP = [code]57133839564458545904789328652610540031895535786011264182548375833179829124845398393126574488675311145377107878746854204162666250198684504466355949195922066574942592095735778929325357290444962472405416790722118445437122269675520000000000000000000000000000000000000 * 19![sup]6[/sup] + 1[/code]? O rly?[/QUOTE]

The '!' was not the factorial-sign!

got it working.

[CODE]sievex(kmin,kmax,b,n,pmax)=a=0;for(k=kmin,kmax,forprime(p=1,pmax,if((k*b!^n+1)%p==0 && (k*b!^n+1)!=p,a=a+1));if(a==0,print("no factor found for"k"*"b"!^"n"+1 in this range"));a=0)[/CODE]

[QUOTE=Karsten]The '!' was not the factorial-sign![/QUOTE]

An integer followed by an exclamation point is taken as a factorial to me.

!

factorial
factorial
combinatorics
n! means the product 1 × 2 × ... × n. 4! = 1 × 2 × 3 × 4 = 24
logical negation
not
propositional logic
The statement !A is true if and only if A is false.

A slash placed through another operator is the same as "!" placed in front.

(The symbol ! is primarily from computer science. It is avoided in mathematical texts, where the notation ¬A is preferred.) !(!A) ⇔ A
x ≠ y ⇔ !(x = y)

according to wikipedia's table.

Yes, yes, I know the connotations of the "!", symbol.

I can cite a few quick examples:

7! = 5040
6 * 6 != 7.

[QUOTE=3.14159;226731]It is a sieve that kicks out bad k-values for numbers of the form k * b![sup]n[/sup] + 1 which are divisible by small primes, where the user is allowed to choose the k-range, the factorial base and the exponent, as well as the largest prime it will sieve up to. The output will be written to a file for testing. The user will have to make the modifications to ensure that the range will be testable, preferably with a decent or advanced text editor that is capable of doing so.

Where k < b![sup]n[/sup], and where k, b, and n are positive integers.[/QUOTE]

:rolleyes:

I've given up trying to explain this. Here's my off-the-top-of-my-head code, with some basic comments:
[code]vk(lim,kmin,kmax,b,n,c=-1)={
my(v=vectorsmall(kmax,j,j>=kmin),start,vv,i=0);
forprime(p=2,lim,
trap(,next, \\ if p can't divide any of the members, just go on the next prime
start=lift(Mod(c,p)/Mod(b,p)^n); \\ start * b^n = c (mod p)
);
start += p*ceil((kmin-start)/p); \\ start >= kmin
forstep(j=start,kmax,p,
v[j]=0 \\ j * b^n = c (mod p)
)
);
vv=vector(sum(i=1,#v,v[i])); \\ vector for the numbers found, rather than 0/1
for(j=1,#v,if(v[j],vv[i++]=j)); \\ fill vector
vv \\ output
};
addhelp(vk, "vk(lim,kmin,kmax,b,n,{c=-1}): Returns those k with kmin <= k <= kmax for which k * b^n - c has no prime divisors up to lim.");[/code]

Test code for k * 19!^6 + 1, k <= 10^4:
[code]vk(1e6,1,10^4,19!,6,-1)[/code]

This takes 200 ms. Comparable trial division takes 1 minute 29 seconds:
[code]B=19!^6;for(k=1,1e4,forprime(p=2,1e6,if((k*B+1)%p==0,break)))[/code]

Impressive job CRG: 200ms for k<10000!

So I've made my first script for WinPFGW, but it uses only trial factoring for small primes and PRP-test the k-values where no factor was found for.

It creates a log-file like this (found primes/PRP are stored in "pfgw.log", too):
[code]
146*19!^6+1: factor 47
147*19!^6+1: factor 31
148*19!^6+1: factor 41
149*19!^6+1: factor 23
150*19!^6+1: is PRP/PRIME.
151*19!^6+1: no factor found.
152*19!^6+1: factor 6269
153*19!^6+1: no factor found.
154*19!^6+1: factor 883
155*19!^6+1: factor 523
156*19!^6+1: factor 65581
157*19!^6+1: factor 113
158*19!^6+1: factor 12583
159*19!^6+1: factor 31707989
[/code]

You can choose min-k, max-k, factorial-n, exponent, +/-1 and factor-bound.
There is no code to prevent using false parameters, but that should not be a problem.

Timing: 1 <= k <= 10000, factor-bound = 100000
The script run 24s (here a AMD XP +2500, 1.8GHz) and found 247 PRP/primes.

wow even with minimal primelimit and trying a near maximal amount before length overflow with my max memory it's working up to 10^7 in under 17 seconds.

[QUOTE=CRGreathouse]I've given up trying to explain this. Here's my off-the-top-of-my-head code, with some basic comments:[/QUOTE]

No need to be rude. I repeatedly said that I was unable to do it on my own.

[QUOTE=science_man_88]wow even with minimal primelimit and trying a near maximal amount before length overflow with my max memory it's working up to 10^7 in under 17 seconds.
[/QUOTE]

It's a nice script, but I still failed in the end. It's a double-edged sword. I can only write basic scripts.

Only 1.5 minutes for 500 million. Nice.

[QUOTE=Karsten]You can choose min-k, max-k, factorial-n, exponent, +/-1 and factor-bound.
There is no code to prevent using false parameters, but that should not be a problem.

Timing: 1 <= k <= 10000, factor-bound = 100000
The script run 24s (here a AMD XP +2500, 1.8GHz) and found 247 PRP/primes.[/QUOTE]

Can you give us the timings for sieving to a certain prime? (Ex: Time it takes to sieve to 10[sup]6[/sup] for a certain range)

Also: Is it directly input to WinPFGW?

[QUOTE=3.14159;226818]Can you give us the timings for sieving to a certain prime? (Ex: Time it takes to sieve to 10[sup]6[/sup] for a certain range)

Also: Is it directly input to WinPFGW?[/QUOTE]

Start WinPFGW.exe and input in the line under "Standard "PFGW compatible" command line" the name of the script-file (in same folder).

That's all.

You can edit the parameters you need in the script.

(Need for sieve to 10^6 and 1<k<10000 about 50 seconds; window minimized!)

Nevermind. I answered my own question. Okay: I can mess with the params as I wish? Cool.

Wait: Does it use the F/trial division switch?

Also: Found 268 * 4200! + 1 (13399 digits). A prime with a prime number of digits. Excellent.

[CODE]plength(n) = forstep(x=(10^n)+1,10^(n+1),[2],if(isprime(x),return(x)))[/CODE]

for some reason this isn't working properly. I can get it to work as planned but then it might mess up someone else if they use it.

realized why now doh.

sm: It looks good to me. This is A003617(n+1), right?

Faster would be
[code]plength(n)=nextprime(10^n)[/code]

[QUOTE=CRGreathouse;226873]sm: It looks good to me.[/QUOTE]

well it isn't! and I know why.

[CODE](15:42) gp > plength(n) = forstep(x=(10^n)+1,10^(n+1),[2],if(isprime(x),return(x)))
%471 = (n)->forstep(x=(10^n)+1,10^(n+1),[2],if(isprime(x),return(x)))
(15:43) gp > plength(3)
%472 = 1009
(15:43) gp >[/CODE]

is the original result

[CODE](15:45) gp > plength(n) = forstep(x=(10^(n-1))+1,10^(n),[2],if(isprime(x),return(x)))
%473 = (n)->forstep(x=(10^(n-1))+1,10^(n),[2],if(isprime(x),return(x)))
(15:45) gp > plength(3)
%474 = 101
(15:45) gp >[/CODE]

see the difference the reason the original was working weird is I didn't account for 1-9 so when i put in plength(3) the first on started at 10^3 which already has 4 digits when I compensate I get the correct value.

I also did this with your knowledge:

[CODE]sumdigits(n) = if(n%9!=0,return(n%9),return(n%9+9))[/CODE]

but you might think it's kinda lame.

[QUOTE]0 [B]1 2 3 4 5 6 7 8 9[/B]
[B]1[/B] [COLOR="Red"]1 2 3 4 5 6 7 8 9[/COLOR]
[B]2[/B] [COLOR="Orange"]2 4 6 8 1 3 5 7 9[/COLOR]
[B]3[/B] [COLOR="Yellow"]3 6 9 3 6 9 3 6 9[/COLOR]
[B]4[/B] [COLOR="green"]4 8 3 7 2 6 1 5 9[/COLOR]
[B]5[/B] [COLOR="Blue"]5 1 6 2 7 3 8 4 9[/COLOR]
[B]6[/B] [COLOR="indigo"]6 3 9 6 3 9 6 3 9[/COLOR]
[B]7[/B] [COLOR="violet"]7 5 3 1 8 6 4 2 9[/COLOR]
[B]8[/B] [COLOR="Magenta"]8 7 6 5 4 3 2 1 9[/COLOR]
[B]9[/B] 9 9 9 9 9 9 9 9 9[/QUOTE]

if i did my math correct this is the sum of digits multiplication table.

[QUOTE=science_man_88;226878]I also did this with your knowledge:

[CODE]sumdigits(n) = if(n%9!=0,return(n%9),return(n%9+9))[/CODE]

but you might think it's kinda lame.[/QUOTE]

This only works for n < 100. For larger numbers, try
[code]dsum(n,b=10)={
my(s=0);
while(n>=b,
s += n%b;
n \= b;
);
s+n
};
addhelp(dsum, "dsum(n,{b=10}): Sum of the base-b digits of n.");[/code]

A more efficient version would work mod 100 or 1000 using a lookup table, and use divrem rather than \ and %. A really efficient version would split the values into two parts of roughly equal length until they became 'small' (say, less than a googol).

[QUOTE=science_man_88;226876]see the difference the reason the original was working weird is I didn't account for 1-9 so when i put in plength(3) the first on started at 10^3 which already has 4 digits when I compensate I get the correct value.[/QUOTE]

Because you didn't explain what the function was supposed to do, you didn't have a help message, and the name wasn't self-explanatory, I couldn't tell that this was an error rather than the intended behavior. :smile:

[QUOTE=3.14159;226844]Also: Found 268 * 4200! + 1 (13399 digits). A prime with a prime number of digits. Excellent.[/QUOTE]

Primes with a prime number of digits... hmm. I wonder in what sense this has positive density. Of course it has natural density 0, as a subset of the primes; I don't know what its upper logarithmic density would be. If that's not positive, is there a finer measure? I wonder.

The set is 'weird' because it has huge gaps, then lots of relatively densely packed elements: no members (strictly) between 10^114 - 153 and 10^126 + 679, a gap of length ~10^126, but then > 3.44218 * 10^123 members in the next 10^126.

I tried it on 2047 it works when did you get an error ?

[QUOTE=science_man_88;226887]I tried it on 2047 it works when did you get an error ?[/QUOTE]

Oh I see: you don't want the sum of the digits but rather the digital sum of the number. In that case your code is fine. (Use addhelp, that way this sort of thing won't happen so much!)

I figured you'd know as it's your advice about it earlier in this thread that I made it with by the way I was able to do:

[CODE](17:21) gp > sumdigits(2^20470000-1)
%550 = 6
(17:21) gp > ##
*** last result computed in 47 ms.[/CODE]

your gave me this:

[CODE](17:21) gp > dsum(2^20470-1)
%549 = 27762
(17:21) gp > ##
*** last result computed in 94 ms.[/CODE]

if we can loop your code until it's a one digit number and speed it up it's okay. but if each step took as long as this one. my code would be 6 times the speed at least.

[QUOTE=CRGreathouse]Primes with a prime number of digits... hmm. I wonder in what sense this has positive density. Of course it has natural density 0, as a subset of the primes; I don't know what its upper logarithmic density would be. If that's not positive, is there a finer measure? I wonder.

The set is 'weird' because it has huge gaps, then lots of relatively densely packed elements: no members (strictly) between 10^114 - 153 and 10^126 + 679, a gap of length ~10^126, but then > 3.44218 * 10^123 members in the next 10^126.[/QUOTE]

Is it as weird as the subset of primes that have a square number of digits? (2, 3, 5, 7, 1009, 1013, 1019, 1021, 1031, 1033, 1039.. 9973, 100000007, ..)

[QUOTE=science_man_88;226889]if we can loop your code until it's a one digit number and speed it up it's okay. but if each step took as long as this one. my code would be 6 times the speed at least.[/QUOTE]

We're computing essentially different things. If all you need is the digital root, then by all means compute that! If you need the sum of the digits, it's a much harder task.

Similarly, factoring an integer is harder than just testing for primality; if you can get away with the latter, do that.


For what it's worth, your code could be modified to be roughly twice as fast. Two possibilities for n > 0 spring to mind:
[code]droot(n)={
my(k=n%9);
if(k,k,9)
};

droot1(n)={
((n-1)%9)+1
};[/code]

I'm not sure which is faster, probably the former. Test them and see!

[QUOTE=3.14159;226890]Is it as weird as the subset of primes that have a square number of digits? (2, 3, 5, 7, 1009, 1013, 1019, 1021, 1031, 1033, 1039.. 9973, 100000007, ..)[/QUOTE]

It's weird in a similar way, yes. Of course things having to do with decimal digits are usually unnatural. :smile:

[QUOTE=CRGreathouse;226891]
I'm not sure which is faster, probably the former. Test them and see![/QUOTE]

by 3 test each the first is ahead by one test as the other 2 match the time for the second. lowest time so far is 16 ms.

[QUOTE=science_man_88;226893]by 3 test each the first is ahead by one test as the other 2 match the time for the second. lowest time so far is 16 ms.[/QUOTE]

The timer isn't really accurate for very small times. Test a lot of numbers to check:
[code]for(i=1,1e5,droot(i))[/code]
or something like that.

Looks to me like droot1 is better, at least for small numbers. For big numbers, let's do
[code]v=vector(10^5,unused,random(10^1000));
for(i=1,#v,droot(v[i]))[/code]

Looks like the winner is still droot1, which is faster than droot, which is faster than sumdigits.

[QUOTE=CRGreathouse]It's weird in a similar way, yes. Of course things having to do with decimal digits are usually unnatural.
[/QUOTE]

Ah well: Random prime of the day:
9845575878878960828055703578149239033062029198602060374047309547974126063126741251854204515517297810571460660189707205143591982239643 (133 digits)

[QUOTE=3.14159;226895]Ah well: Random prime of the day:
[code]
9845575878878960828055703578149239033062029198602060374047309547974126063126741251854204515517297810571460660189707205143591982239643 (133 digits)[/code][/QUOTE]

...and
[code]
98455758788789608280557035781492390330620291986020603740473095479741260631267412518542045155172978105714606601897072051435919822396433
[/code]
(134 digits), too!

[QUOTE=Karsten]...and

[code]98455758788789608280557035781492390330620291986020603740473095479741260631267412518542045155172978105714606601897072051435919822396433[/code]

(134 digits), too!
[/QUOTE]

[code]160074594123956794219458402545764711419701756584167638068329365693034831201396939643913504262352130745852290381939420288140978693341478528259387657797580376902986395224408345354821071590949525350364018649578862749297255505837416439606095853836355923951566526457673735312339184330365326723420347366112675200713878660771081227364603472040699601553789595479339960277304318254798850917751885704667910151149091532150905626179048799036700679933820859845365985390911277048117168671365128451668210581621327277787092273718680191991994664786587234396792965881059281805213980760310898060678555989602303590308736441654087832242397146458759717551025723008092342436218853265799505252096172021953238850043510060599858403400384570222716416699106818277439891613536333361593153540409685328569171314344686601795619300473375554423332099932225456784121921559118525805374656294079430382561204188488274552102452626811652972822901469815861523877416162913528086735710854833100692524526048253371898482177626664655396682221293652089791416114579511127188674154339902564205267260071373223290668469512503977940823503507668484572724611665887674632905471958304064277923482051286050072649[/code]

1155 digits! Ha!

[code]160074594123956794219458402545764711419701756584167638068329365693034831201396939643913504262352130745852290381939420288140978693341478528259387657797580376902986395224408345354821071590949525350364018649578862749297255505837416439606095853836355923951566526457673735312339184330365326723420347366112675200713878660771081227364603472040699601553789595479339960277304318254798850917751885704667910151149091532150905626179048799036700679933820859845365985390911277048117168671365128451668210581621327277787092273718680191991994664786587234396792965881059281805213980760310898060678555989602303590308736441654087832242397146458759717551025723008092342436218853265799505252096172021953238850043510060599858403400384570222716416699106818277439891613536333361593153540409685328569171314344686601795619300473375554423332099932225456784121921559118525805374656294079430382561204188488274552102452626811652972822901469815861523877416162913528086735710854833100692524526048253371898482177626664655396682221293652089791416114579511127188674154339902564205267260071373223290668469512503977940823503507668484572724611665887674632905471958304064277923482051286050078929[/code]
1155 digits. :lol:

CRG is there a way to create a Vec(Vec(input())) such that each element is a string if so we can search it word by word instead of letter by letter for something to help you and simply replace one for another.

[QUOTE=science_man_88;226909]CRG is there a way to create a Vec(Vec(input())) such that each element is a string if so we can search it word by word instead of letter by letter for something to help you and simply replace one for another.[/QUOTE]

I don't know what you're looking for, exactly, but that sounds inefficient.

What about doing this instead:
[code]substring(string_to_search, string_to_find)={
my(s1=Vec(string_to_search),s2=Vec(string_to_find),good);
for(i=0,#s1-#s2,
good=1;
for(j=1,#s2, if(s2[j] != s1[j+i], good=0;break));
if(good, return(1))
);
0
};[/code]

[code]> substring("good day", "day")
%1 = 1
> substring("good day", "night")
%2 = 0[/code]

Of course there are better ways to do this, but I wouldn't worry about that until performance is an issue. At that point you can either make it more efficient or you can use a different system entirely.

see instead of returning 1 or 0 I was going to replace the string I found with what Pari can understand. but I'll try it your way for now.

is there a way to replace what we find with another string ? if so maybe say if this is found return("stringx") into the string instead.

[QUOTE=science_man_88;226919]is there a way to replace what we find with another string ? if so maybe say if this is found return("stringx") into the string instead.[/QUOTE]

Sure, make whatever changes you like.

I think we can use what we did for shifting to transfer from position j+ length(s2) = to a new Vec then delete the s2 from where it was found concat(newstring) on the new end of it then concat (the rest) back on to the string after adjusting j for the new length of the string so it only checks what it hasn't before. then repeat for other strings. then we can use these names placed in to print the function by name and if we can get it to construct it in the correct manner from there it'll be done except how to read into the s1.

Why so complicated?

Because it is PARI -> program to calculate mathematic formulas and values or show functions, but not to manipulate texts like search/replace strings and other string-functions!

What you want to create is hard enough to be successful, but using PARI for that is even impossible!

Taking another tool/language would be better! Trust me.

[QUOTE=kar_bon;226928]Why so complicated?

Because it is PARI -> program to calculate mathematic formulas and values or show functions, but not to manipulate texts like search/replace strings and other string-functions!

What you want to create is hard enough to be successful, but using PARI for that is even impossible!

Taking another tool/language would be better! Trust me.[/QUOTE]

our last conversation on this was post 285 kar_bon up to you if you want the task of writing it in another language meanwhile I'm learning about how to make Pari do work for CRG lol.

[QUOTE=science_man_88;226929]our last conversation on this was post 285 kar_bon up to you if you want the task of writing it in another language meanwhile I'm learning about how to make Pari do work for CRG lol.[/QUOTE]

Ok, I wish you luck for that task! And patience!

Worthless analogy of the day:

sm is building what looks like a shed, pounding in nails with a coconut. kar_bon says, "Instead of using that coconut, why not use this hammer?". CRGreathouse says, "Good start, but making a rocketship is hard, even with a hammer." sm says, "You don't get it, I'm building a rocketship with this coconut because I want to get better with coconuts (and because I want a rocketship).".

[QUOTE=CRGreathouse]Worthless analogy of the day:

sm is building what looks like a shed, pounding in nails with a coconut. kar_bon says, "Instead of using that coconut, why not use this hammer?". CRGreathouse says, "Good start, but making a rocketship is hard, even with a hammer." sm says, "You don't get it, I'm building a rocketship with this coconut because I want to get better with coconuts (and because I want a rocketship).".
[/QUOTE]

What an accurate analogy.

fine then WTF do you get me interested ? if you don't want me to do anything with life don't give me life, if you don't want me to do within Pari then why did you bother to get me into it in the first place and WTF is the forum here if you constantly seem to not care what others think.

[QUOTE=science_man_88;226940]fine then WTF do you get me interested ? if you don't want me to do anything with life don't give me life, if you don't want me to do within Pari then why did you bother to get me into it in the first place and WTF is the forum here if you constantly seem to not care what others think.[/QUOTE]

I think it's great that you (and pi) are learning Pari, and I don't mind if you aren't able to finish your project. I've certainly started projects that I wasn't able to finish, but were worthwhile.

I see two reasonable approaches here. There's the KB plan: learn C or Perl and write your parser. The parser won't be perfect, but it will be able to transform some queries into valid Pari code, and that's pretty cool.

The other approach is the SM plan: keep writing it in Pari. The parser won't be as easy to write, but you'll learn a whole lot more about Pari than you would with the other method.

Of course you could also give up, but in that case you have neither parser nor Pari skills. :smile:

lets put it this way I know more syntax in Pari then i do in any other language I've learned before , so me doing Perl or C is likely going to never work. you seem to be the only one I've talked to have an interest other than me in the PARI method and I sent a few PM's to check for someone who knows perl but I haven't found anyone.

[QUOTE=science_man_88]fine then WTF do you get me interested ? if you don't want me to do anything with life don't give me life, if you don't want me to do within Pari then why did you bother to get me into it in the first place and WTF is the forum here if you constantly seem to not care what others think.
[/QUOTE]

Stop raging, this is merely an Internet message board for math enthusiasts. Follow this: "As soon as I exit this forum, you no longer exist."

[QUOTE=3.14159;226948]Stop raging, this is merely an Internet message board for math enthusiasts. Follow this: "As soon as I exit this forum, you no longer exist."[/QUOTE]

I rarely have the forum closed during the day even when I'm away I have it open most days.

[QUOTE=science_man_88]I rarely have the forum closed during the day even when I'm away I have it open most days.
[/QUOTE]

[B]Last off-topic post[/B]: "As soon as I stop viewing this forum, you no longer exist." ?

[QUOTE=3.14159;226948]Stop raging, this is merely an Internet message board for math enthusiasts. Follow this: "As soon as I exit this forum, you no longer exist."[/QUOTE]

[QUOTE=science_man_88;226949]I rarely have the forum closed during the day even when I'm away I have it open most days.[/QUOTE]

That's the best pair of snarky posts I've seen all week. :grin:

Now: Where were we?

[quote=CRGreathouse]That's the best pair of snarky posts I've seen all week. :grin:[/quote]

Well, only hopeless failures, internet addicts, and/or gullible morons believe someone actually cares for them over the Internet. It's simply the truth.

[QUOTE=3.14159;226953]Now: Where were we?



Well, only hopeless failures, internet addicts, and/or gullible morons believe someone actually cares for them over the Internet. It's simply the truth.[/QUOTE]

you've defined me and my life dead on. happy days.

[URL="http://www.numberempire.com/texequationeditor/equationeditor.php"]Here's[/URL] a good place if you want to fiddle around with TeX. (LaTeX, to be exact. I don't know the differences between the two.)

[CODE]forstep(x=#v-1,1,[-1],for(i=1,#v-1,[COLOR="Red"]v[i]=v[i+1]-v[i][/COLOR]);v=vector(#v-1,i,v[i]))[/CODE]

came with this today as I lost what i had last night looks like my computer got restarted by my mom as she doesn't know how to turn it from full sleep to awake and with all the data intact. I have made it past this again but I'm having trouble trying to make it figure out the equation lol.

Well, I don't know what the code is supposed to do, but I'll have a look.

[CODE]forstep(x=#v-1,1,-1,
for(i=1,#v-1,
v[i]=v[i+1]-v[i]
);
v=vector(#v-1,i,v[i])
)[/CODE]

First problem: you're using values before you define them. If, mathematically, [TEX]v_i=v_{i+1}+v_i[/TEX], then you can rearrange this as [TEX]v_{i+1}=0[/TEX] and do
[CODE]v[1]=foo;
for(i=2,#v,
v[i]=0
);
v[/CODE]
on the inner loop, where foo is the initial value (I don't know what this would be).

The second problem is that the outer loop doesn't do anything with the values of v that you calculate, and the inner loop doesn't use the x.

This makes your code essentially
[code]nop(x)={

};[/code]

The third problem is that you're using v without defining it, twice: once in the outer loop, for the starting value of x, and then in the inner. So actually the code is more like
[code]foo(x)={
if(type(v) == "t_VEC" || type(v) == "t_COL",
for(i=1,#v-1,v[i]=v[i+1]-v[i]);
v=vector(#v-1,i,v[i])[COLOR="Red"];[/COLOR]
,
error("_[_]: not a vector.")
)
};[/code]
where the red semicolon emphasizes that this is *not* output (it's a side-effect).

okay I'll add a V but I wanted it user defined before hand. it calculates sequence of difference then.

maybe I didn't give it a name because i thought you may misinterpret it as you did my last attempt at anything.

[QUOTE=science_man_88;227026]okay I'll add a V but I wanted it user defined before hand. it calculates sequence of difference then.[/QUOTE]

That's a good idea, but in this case make it a function argument:
[code]foo(v,x)={
\\ stuff...
};[/code]

[QUOTE=science_man_88;227026]maybe I didn't give it a name because i thought you may misinterpret it as you did my last attempt at anything.[/QUOTE]

It's a funny thing, as you get better and better with Pari the most difficult part is now usually just explaining what you want rather than coding it. :smile:

[QUOTE=CRGreathouse;227036]That's a good idea, but in this case make it a function argument:
[code]foo(v,x)={
\\ stuff...
};[/code]



It's a funny thing, as you get better and better with Pari the most difficult part is now usually just explaining what you want rather than coding it. :smile:[/QUOTE]

hilarious! not. I think for the finding the equation i can store the value i get at the level as coefficients in a Vec and subtract it from a vector like the original v and then do the loop again so this may take another loop yet got a faster way ?

[QUOTE=science_man_88;227037]I think for the finding the equation i can store the value i get at the level as coefficients in a Vec and subtract it from a vector like the original v and then do the loop again so this may take another loop yet got a faster way ?[/QUOTE]

I don't know what you're doing, so I can't speak to whether this is a correct or an efficient way to do what you want. But Pari can do this, certainly.

sequence of difference like A134458 is for A000043 but as deep as it can to find an equation for the sequence.

[QUOTE=science_man_88;227040]sequence of difference like A134458 is for A000043 but as deep as it can to find an equation for the sequence.[/QUOTE]

Ah! Got it.

[code]diff(v)={
vector(#v-1,i,v[i+1]-v[i])
};
addhelp(diff, "diff(v): First differences of the vector v.");[/code]

Sample:
[code]> diff([2,3,5,7,11,13,17,19,23])
%1 = [1, 2, 2, 4, 2, 4, 2, 4][/code]

[QUOTE=science_man_88;227040]to find an equation for the sequence.[/QUOTE]

I have some code for something like that. It's ugly, and it doesn't always work, but try this:

[code]findrec(v:vec, verbose:bool=1)={
my(c,d = (#v - 1) >> 1, firstNonzero = 0);
if (#v < 3,
error("findrec: Not enough information to determine if the sequence is a recurrence relation: matrix is underdetermined. Give more terms and try again.")
);
forstep(i=d,1,-1,
if(v[i] != 0, firstNonzero = i)
);
if(firstNonzero == 0,
if(verbose,
print1("Recurrence relation is a(n) = 0.");
);
return([0]~);
);
for (i=firstNonzero,d,
c = findrecd(v,i,verbose);
if(c, return(c))
);
if(verbose,print("Cannot be described by a homogeneous linear recurrence relation with "d" or fewer coefficients."));
0
};
addhelp(findrec, "findrec(v, {verbose=1}): Tries to find a homogeneous linear recurrence relation with constant coefficients to fit the vector v.");


findrecd(v:vec, d:int, verbose:bool=1)={
my(M,c);
if (#v < 2*d,
error("findrec: Not enough information to determine if the sequence is a "d"-coefficient recurrence relation; matrix is underdetermined. Give more terms or reduce d and try again.")
);
M = matrix(d,d,x,y,v[d+x-y]);
\\print(M" * c = "vector(d,i,v[d+i])~);
c = matsolve(M,vector(d,i,v[d+i])~);
for(n=2*d+1,#v,
if(v[n] != sum(i=1,d,v[n-i] * c[i]),return(0))
);
if(verbose,
my(init=1,s);
print1("Recurrence relation is a(n) = ");
for(i=1,#c,
if(c[i] == 0, next);
if(init,
s = initial(c[i], Str("a(n-", i, ")"));
init = 0
,
s = Str(s, medial(c[i], Str("a(n-", i, ")")))
)
);
print(s".");
if((vecmax(c) == 1 && vecmin(c) == 0 && vecsum(c) == 1) || c == [1]~,
print("Sequence has period "#c".");
,
my(g=0);
for(i=1,#c,
if(c[i] != 0, g = gcd(g, i))
);
if (g > 1,
my(gvec = vector(#c/g, i, c[i*g]),s,init=1);
for(i=1,#gvec,
if(gvec[i] == 0, next);
if(init,
s = initial(gvec[i], Str("a(n - ", i, ")"));
init = 0
,
s = Str(s, medial(gvec[i], Str("a(n - ", i, ")")))
)
);
print("Can be though of as "g" interlocking sequences, each of the form a(n) = "s".")
)
);
print((#v-2*d)" d.f.")
);
c
};
addhelp(findrecd, "findrecd(v, d, {verbose=1}): Helper function for findrec. Tries to find a d-coefficient homogeneous linear recurrence relation with constant coefficients to fit the vector v.");

initial(n:int, s:str)={
if (n == 0, return(""));
if (n == -1, return(Str("-", s)));
if (n == 1, return(s));
Str(n, s)
};


medial(n:int, s:str)={
if (n == 0, return(""));
if (n == -1, return(Str(" - ", s)));
if (n == 1, return(Str(" + ", s)));
Str(if (n < 0, " - ", " + "), abs(n), s)
};[/code]

Example:
[code]> findrec([1,1,2,3,5,8,13])
Recurrence relation is a(n) = a(n-1) + a(n-2).
3 d.f.
%1 = [1, 1]~[/code]

too bad it doesn't seem to recognise [1,3,7] as anything.

[QUOTE=science_man_88;227048]too bad it doesn't seem to recognise [1,3,7] as anything.[/QUOTE]

That's too short for it to be able to do much analysis. It could construct a recurrence relation for it (with two coefficients), but it's programmed to not do that, since there would be no more terms to check that the relation is correct and not just spurious.

For example, it fits the relation
6a(n-1) - 11a(n-2)
and the relation
5a(n-1) - 8a(n-2)
etc.

good now I can check Mersenne exponents for things to eliminate.

[QUOTE=science_man_88;227051]good now I can check Mersenne exponents for things to eliminate.[/QUOTE]

On my main computer I can just type
[code]findrec(bfile(43))[/code]
where 43 means A000043. :smile:

[QUOTE=CRGreathouse;227052]On my main computer I can just type
[code]findrec(bfile(43))[/code]
where 43 means A000043. :smile:[/QUOTE]

I was trying to find patterns in what to eliminate but cool. :lol:

though I'll have to look at this closer [url]http://www.research.att.com/~njas/sequences/A055010[/url]

every prime above 5 in this sequence seems fit to eliminate.

and I tried the same basic thing (which gave me the same equation by your code) starting on 29 as well but I haven't checked that one either.

[QUOTE=science_man_88]I was trying to find patterns in what to eliminate but cool.
[/QUOTE]

Mersenne numbers do not have any patterns, they're randomly distributed and have no covering set of divisors (Besides having Proth numbers as factors)