mersenneforum.org

mersenneforum.org (https://www.mersenneforum.org/index.php)
-   YAFU (https://www.mersenneforum.org/forumdisplay.php?f=96)
-   -   Yafu (https://www.mersenneforum.org/showthread.php?t=10871)

bsquared 2011-12-01 14:22

[QUOTE=LaurV;280641]now another one:

try this:

yafu factor(460707809)
[/QUOTE]

Aha! This factorization is via fermat, and I had apparently forgot to include code to log factors found via that method. It will be an easy fix.

Give yourself a merit badge in beta testing.

bsquared 2011-12-01 14:50

[QUOTE=jasonp;279707]
PS: Ben, do you compute the number of digits in the decimal representation of a number by taking the number of bits and scaling? Msieve used to do that too, but when enough people complained that this was occasionally off by a digit I decided to just print the damn number to a string and measure its size :)[/QUOTE]

I used to repeatedly divide by 10 and count the number of times I divided (silly, I know, but it worked). Now, I just call mpz_sizeinbase(num,10). So the problem is in mpz_sizeinbase???

Maybe sizeinbase uses bits and scaling.

bsquared 2011-12-01 15:36

v. 1.29.2 available
 
Includes fixes for the most recently reported bugs.

Also removed zlib support, so data files will be as before.

[URL="https://sourceforge.net/projects/yafu/"]Link[/URL].

axn 2011-12-01 15:39

[QUOTE=bsquared;280668]I used to repeatedly divide by 10 and count the number of times I divided (silly, I know, but it worked).[/quote]
Really? Not 10^9 (32-bit) or 10^19 (64-bit)?

[QUOTE=bsquared;280668]Now, I just call mpz_sizeinbase(num,10). So the problem is in mpz_sizeinbase???

Maybe sizeinbase uses bits and scaling.[/QUOTE]

[url]http://www.delorie.com/gnu/docs/gmp/gmp_42.html[/url] says [QUOTE]Return the size of op measured in number of digits in base base. The base may vary from 2 to 36. The sign of op is ignored, just the absolute value is used. The result will be exact [B]or 1 too big[/B][/QUOTE] :shock:

bsquared 2011-12-01 16:17

[QUOTE=axn;280674] Really? Not 10^9 (32-bit) or 10^19 (64-bit)?
[/QUOTE]

Yup, really. Just one of those things I did just to get done.

But, since it's fun and I'm waiting for something else to compile anyway, how's this:

[CODE]

[FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]int[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] ndigits_1(fp_digit n)[/SIZE][/FONT]
[SIZE=2][FONT=Consolas]{[/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff] int[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] i=0;[/SIZE][/FONT]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff] while[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] (n != 0)[/SIZE][/FONT]
[SIZE=2][FONT=Consolas] {[/FONT][/SIZE]
[SIZE=2][FONT=Consolas] n /= 10;[/FONT][/SIZE]
[SIZE=2][FONT=Consolas] i++;[/FONT][/SIZE]
[SIZE=2][FONT=Consolas] }[/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff] if[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] (i==0)[/SIZE][/FONT]
[SIZE=2][FONT=Consolas] i++;[/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff] return[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] i;[/SIZE][/FONT]
[SIZE=2][FONT=Consolas]}[/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff]
[SIZE=2][FONT=Consolas][COLOR=#0000ff]int[/COLOR][/FONT][/SIZE][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] ndigits(z *n)[/SIZE][/FONT]
[SIZE=2][FONT=Consolas]{[/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff] int[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] i=0;[/SIZE][/FONT]
[SIZE=2][FONT=Consolas] z nn,tmp;[/FONT][/SIZE]
fp_digit r;
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#008000][FONT=Consolas][SIZE=2][COLOR=#008000][FONT=Consolas][SIZE=2][COLOR=#008000] //can get within one digit using zBits and logs, which would[/COLOR][/SIZE][/FONT]
[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#008000][FONT=Consolas][SIZE=2][COLOR=#008000][FONT=Consolas][SIZE=2][COLOR=#008000] //be tons faster, but which is sometimes wrong. this is slower [/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT]
[FONT=Consolas][SIZE=2][COLOR=#008000][FONT=Consolas][SIZE=2][COLOR=#008000][FONT=Consolas][SIZE=2][COLOR=#008000] //but exact.[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT]
[FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] zInit(&nn);[/SIZE][/FONT]
[SIZE=2][FONT=Consolas] zInit(&tmp);[/FONT][/SIZE]
[SIZE=2][FONT=Consolas] zCopy(n,&tmp);[/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff] while[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] (tmp.size > 1)[/SIZE][/FONT]
[SIZE=2][FONT=Consolas] {[/FONT][/SIZE]
[SIZE=2][FONT=Consolas] zCopy(&tmp,&nn);[/FONT][/SIZE]
[SIZE=2][FONT=Consolas] r = zShortDiv(&nn,MAX_DEC_WORD,&tmp);[/FONT][/SIZE]
[SIZE=2][FONT=Consolas] i += DEC_DIGIT_PER_WORD;[/FONT][/SIZE]
[SIZE=2][FONT=Consolas] }[/FONT][/SIZE]
[SIZE=2][FONT=Consolas] i += ndigits_1(tmp.val[0]);[/FONT][/SIZE]
[SIZE=2][FONT=Consolas] zFree(&nn);[/FONT][/SIZE]
[SIZE=2][FONT=Consolas] zFree(&tmp);[/FONT][/SIZE]
[/SIZE][/FONT][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff][FONT=Consolas][SIZE=2][COLOR=#0000ff] return[/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][/COLOR][/SIZE][/FONT][FONT=Consolas][SIZE=2][FONT=Consolas][SIZE=2] i;[/SIZE][/FONT]
[SIZE=2][FONT=Consolas]}[/FONT][/SIZE]
[/SIZE][/FONT][/CODE]

[QUOTE=axn;280674]
:shock: [/QUOTE]

Ditto.

Andi47 2011-12-01 20:45

#digits = 1+floor(log[sub]10[/sub](n)) would not work?

bsquared 2011-12-01 20:55

[QUOTE=Andi47;280698]#digits = 1+floor(log[sub]10[/sub](n)) would not work?[/QUOTE]

No, because I don't have any way to directly evaluate log10(n) (keep in mind n could be a thousand digits long). So I would compute log2(n) * log10(2) because log2 is easy in binary. Using log2 is what causes the inaccuracy - 99 and 100 both have 7 bits, for instance (and thus would both evaluate to 3 using your formula).

What I could do, I suppose, is use the above trick to estimate ndigits, then compute 10^(ndigits) to make sure it's smaller than the input. That would be faster than iteratively dividing by 10, but we'd be splitting hairs.

EdH 2011-12-02 03:39

I'm totally out of my knowledge base, so bear that in mind, please. Could you use something as simple as:
[code]
if (n < 10^x)
d--;
[/code]??

OK, I'm going away now, for a little while.:smile:

bsquared 2011-12-02 05:27

[QUOTE=EdH;280726]I'm totally out of my knowledge base, so bear that in mind, please. Could you use something as simple as:
[code]
if (n < 10^x)
d--;
[/code]??

OK, I'm going away now, for a little while.:smile:[/QUOTE]

Yes, something similar, although you'd need to say how you get x and d. What I've now got implemented is essentially the following:

[CODE]
guess = ceil(log2(n) * log10(2));
temp = 10 ^ (guess - 1)
if (temp > n)
guess--;
return guess;
[/CODE]

EdH 2011-12-02 17:39

Your version seems to be the same as mine, programmatically, with your guess = to my d.

I should have explained further, but actually figured x would be the current value for digits minus 1, digits being d. Using the present value for digits as d, I suppose the better way to code it would be:
[code]
if(n < 10^(d-1))
d--;
[/code]n is the value of the number
d is the currently calculated digits (that are giving you an occasional +1)

So, for your example of 99, you get 3 digits (erroneously).

So you compare 99 to 10^(3-1)=100 and since it is less, decrement d.
d becomes 2 digits whenever it's 1 too many.

Since you are only ever off by one and the error is always positive, the above should be all that is necessary to fix it. But, the question is why isn't yours working?

Of course, I don't know the intricacies of handling the really big numbers...

jasonp 2011-12-02 23:12

Ben, you could also use a function like mp_log in Msieve, that limits the evaluation of the input to its top 2 or 3 words and then multiplies by the number of bits below that.


All times are UTC. The time now is 22:51.

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.